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1. 


INTRODUCTION 


Classical  initial  value  problems  of  systems  of  first-order  ordinary  differential 
equations  are  designed  to  solve 


y'=f(t,y);  y(0)  =  y 


(i.i) 


Conventional  methods  such  as  Runge-Kutta  and  linear  multistep  (LMS)  methods 
have  been  very  well  developed  (Henrici  [16])  and  have  been  demonstrated  effective 


for  solving  (1. 1)  when 


df 


ay 


is  small.  Frequently,  for 


af 


*y 


large,  prohibi¬ 


tively  small  step  size  values  (h)  are  required  for  accuracy;  then,  conventional 
methods  seem  impractical.  To  overcome  this  difficulty,  we  choose  to  write  (1. 1) 
in  the  following  way: 


y'  -  Ay  +  g(t»y);  y(°)~y0*  (1*1)' 

where  f(t.y)  may  be  written  as  Ay +  g(t,y)  and  A  is  either  a  constant  matrix  or 
a  function  of  t.  In  the  case  where  Re{x(A)}  <  0  and  X(A),  the  eigenvalues  of  A, 
differ  greatly  in  magnitude  and  g(t,y)  is  a  slowly  varying  function  in  t,  equation 
(1. 1)’  is  called  a  "stiff"  equation.  These  stiff  equations  frequently  occur  in  the 
applications  to  chemical  kinetics,  reactor  calculations,  missile  guidance,  etc. 

The  search  for  effective  schemes  to  solve  stiff  equations  began  over  two 
decades  ago  and  still  goes  on.  Curtiss  &  Hirshfelder  [9]  encountered  the  stiff 
phenomenon  in  their  study  of  chemical  kinetics  and  proposed  low-order  multistep 
formulas  to  integrate  scalar  stiff  equations.  Cohen  [8],  in  solving  reactor  kinetics 
equations,  presented  a  generalization  of  Runge-Kutta  methods.  Certaine  [7] 
demonstrated  that  if  conventional  schemes,  such  as  trapezoidal  rule,  were  used  to 
solve  (1. 1)’,  then  two  problems  were  encountered  —  step  size  and  accuracy. 
Certaine  then  proposed  a  method  to  handle  scalar  stiff  equations  that  have  short 
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time  constants.  Dahlquist  [10]  discussed  a  general  treatment  of  the  stability  of 
linear  multistep  methods  and  investigated  the  special  stability  problem  in  connec¬ 
tion  with  stiff  equations.  Dahlquist  introduced  an  important  concept,  A-stability, 
and  proved  that  there  do  not  exist  A-stable  methods  among  linear  multistep  methods 
of  order  higher  than  2.  Widlund  [32]  and  Gear  [  14]  relaxed  the  A-stability  con¬ 
cept  in  an  attempt  to  create  higher  order  multistep  formulas  suitable  for  stiff  equa¬ 
tions.  Widlund  [32]  defined  A(a) -stability  and  showed  that  there  exist  K-step 
methods  of  order  K  which  are  A(ar)-stable  for  any  a  <  w/2  and  K  <  4.  Gear  [14] 
weakened  the  A-stability  concept,  defined  stiff- stability,  and  derived  stiffly  stable 
methods  of  order  <6.  Gear  [14]  designed  a  computer  program  to  perform  such 
calculations.  Stiffly  stable  methods  of  orders  7  and  8  have  been  found  by  Dill  (1969) 
and  of  order  up  to  11  by  Jain  (1970),  but  no  tests  have  been  made  on  their  algo¬ 
rithms.  Lawson  [24]  generalized  the  Runge-Kutta  method,  Norsett  [28]  general¬ 
ized  the  Adams -Bashforth  methods,  and  Bjurel  [3]  modified  linear  multistep 
methods.  Contributions  also  have  been  made  by  Miranker  [27]  and  Guderley 
et  al.  [15],  who  designed  stiff  methods  to  solve  special  types  of  stiff  equations. 

The  theory  for  stiff  systems  lacks  a  cohesiveness  that  this  thesis  attempts 
to  achieve  by  providing: 

(1)  A  complete  formulation  of  nonlinear  multistep  (NLMS)  methods,  which 
is  demonstrated  to  be  a  generalization  of  linear  multistep  methods  both 
in  technique  and  in  theory. 

(2)  A  full  development  and  proof  of  the  theory  of  NLMS  methods  with 
regard  to  stability,  consistency,  and  convergence. 
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(3)  A  proof  that  NLMS  methods  accommodate  A-stability  in  the  sense  of 
Dahlquist  [10], 

(4)  A  study  of  the  effect  of  the  error  function  Cp+^  (as  defined  in  section 
4,  Computational  Considerations)  by  means  of  a  perturbation  of  the 
characteristic  roots.  The  study  shows  that  LMS  methods  of  order  p  that 
possess  the  smallest  error  C  p+i  are  only  weakly  stable.  However, 
it  will  be  indicated  that  there  always  exists  a  NLMS  family  possessing 
the  smallest  error  C  ... 

p+1 

(5)  Extensive  tests  of  NLMS  methods  applied  to  a  set  of  selected  scalars 
and  systems  of  stiff  equations.  Results  are  compared  with  Adams' 
methods  [16],  Gear's  program  [14],  Seinfeld's  paper  [29]  and  Ehle's 
research  report  [ll] ,  and  it  is  shown  that  NLMS  has  definite  advan¬ 
tages  over  the  above  techniques. 

(6)  A  section  of  conclusions  and  a  summary  of  remaining  problems  with 
some  suggested  solutions. 

(7)  A  listing  of  computer  programs  used  to  implement  the  NLMS  methods. 


2. 


PRELIMINARY  CONSIDERATIONS 


In  this  section,  we  define  the  problems  under  consideration  and  state  the 
theory  in  relation  to  the  existence  and  uniqueness  of  the  solutions  of  approximating 
difference  equations;  the  proofs  for  the  existence  and  uniqueness  theorem  can  be 
found  in  references  [16]  and  [18].  The  starting  procedure  involves  the  use  of 
initial  data;  the  solution  of  the  difference  equations  depends  continuously  on  the 
initial  data.  The  order  of  the  multistep  methods  will  be  defined  as  they  are  dis¬ 
cussed  and  developed  in  section  3,  Theory. 

2. 1.  Problems  Considered 

In  this  paper,  we  consider  the  initial  value  problems  of  a  system  of  first- 
order  ordinary  differential  equations  of  the  form- 

y'  =  Ay  +  fl(t.y) 

=  f(t,y)  ;  y(0)=yo  <2.i) 

in  the  region  R,  defined  byO  =  a<t<b<«°;  ||y||<».  A  is  either  a  constant 
matrix  or  a  function  of  t;  consequently,  a  portion  of  the  theory  will  be  restricted 
to  the  important  case,  where  the  differential  equations  are  stiff,  i.  e. ,  Re{\(A)}  <  0. 
Among  the  numerical  test  examples,  A  is  chosen  to  be  a  constant  matrix  and 
Re{x(A)}  <  0*  The  function  g(t,y)  C  C^1  (p  >  0)  satisfies  the  Lipschitz  condition, 

||fl(t.y*)  -  g(t,y)||  <  L  ||y*  -  y|| .  (2. 2) 

For  the  most  interesting  applications  (those  restricting  the  step  size  to  conven¬ 
tional  methods)  p(A)  »  L,  where  p(A)  is  the  spectral  radius  of  A. 
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2.  2  Existence  and  Uniqueness  Theorem 

We  assume  that  our  initial  value  problems  satisfy  the  conditions  required 
by  the  existence  and  uniqueness  theorem.  We  state  the  existence  and  uniqueness 
theorem  expressed  with  respect  to  f(t,y)  as  follows: 

Theorem  2.  2:  (Existence  and  Uniqueness  Theorem) 

We  assume  that  f(t,y)  satisfies  the  following  two  conditions: 

(1)  f(t,y)  is  continuous  in  R,  where  R  is  the  region 
0  =  a<t<b<«c,|[y||  <«o. 

(2)  3  a  Lipschitz  constant  L*  9  for  arbitrary  t£[a,b]  and  any  two 
vectors  y  and  y*,  the  following  condition  is  satisfied: 

||l(t,y*)-f(t,y)||<L*||y*-y||.  (2.3) 

Then,  for  any  given  initial  vector  yQ ,  3  one  and  only  one  y(t)  ? 

(1)  y(t)  is  continuous  and  continuously  differentiable  for  t£[a,b] 

(2)  y'(t)  =  f  (t,y),  t  C  [a, b] 

(3)  y(a)  =  y0- 

To  the  differential  equations  (i.  1),  we  adjoin  appropriate  relations,  called 
initial  conditions,  that  serve  to  define  a  "meaningful  problem.  "  If  the  solution  of 
(1. 1)  satisfies  appropriate  initial  conditions  of  smoothness,  the  problem  (1. 1)  is 
termed  well  posed  in  the  sense  of  Hadamard  (Isaacson  [20]),  and  the  problem  has 
a  bounded,  unique  solution.  Hochstadt  [18]  proved  that  the  solution  of  the  differ¬ 
ential  equation  depends  continuously  on  the  initial  data.  This  then  fulfills  the 
Hadamard  well-posed  statement.  It  ought  to  be  pointed  out  that  even  though 
Hadamard' s  well-posed  criterion  is  fulfilled,  conventional  techniques  fail  where 
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there  is  a  very  large  Lipschitz  constant.  This  is  the  area  where  we  need  to  con¬ 
sider  nonlinear  multistep  methods. 


2.3.  Norms 

In  a  finite  n-dimensional  vector  space,  we  define  the  p-norm  of  a  vector 

to  be 


n  xl/p 

r*  i  ~  ip\ 


w 

jxii-E  i*,i 


#xL"^2'„»xHp'nfIlxjl- 

For  a  matrix  A  of  order  n,  we  denote 

A  (A)  as  the  eigenvalues  of  A  and 
pCA)  as  the  spectral  radius  of  A  . 
The  different  norms  of  A  take  the  following  definitions: 

llANs  =[»<A*A))l/2 


n 

max  laJ 

n 

max  la  I 

i<i<n j=i 1  y 


^  2  \1/2 


M££-i 
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Let  ||  •  ||^  indicate  one  norm  and  let  ||  •  ||g  indicate  another  norm.  Then 
II  *  II a  and  II  *  II q  are  said  to  be  equivalent  if  3  two  positive  numbers  a  and  b  * 


We  know  that  in  a  finite-dimensional  space,  all  norms  are  equivalent;  therefore,  the 
norms  used  in  this  paper  do  not  refer  to  any  specific  norm.  However,  in  the  Pade 
approximation,  we  use  the  column  norm  ||  •  ||^  ;  in  the  Ehle's  test  examples,  we  use 
||  •  || 2 ;  and  in  the  other  test  examples,  we  use  II  *11^.  In  the  lemmas  and  theorems 
where  the  norm  does  not  have  a  subscript  we  mean  any  norm. 


X 


3. 


THEORY 


This  section  develops  the  full  theory  that  provides  the  basis  for  NLMS 
methods,  including  the  development  of  the  important  notions  of  Stability  and  Con¬ 
sistency.  The  formulation  of  NLMS  methods,  which  assures  consistency  and  is 
closely  connected  with  the  theory,  is  also  described  in  this  section. 

NLMS  methods  are  demonstrated  to  be  a  generalization  of  LMS  methods. 

The  general  formula  is  described  by  (3.5),  and  the  formulation  is  expressed  by 
(3. 36).  Both  explicit  and  implicit  schemes  are  given  in  matrix  expressions  up  to 
an  order  of  3.  The  application  of  NLMS  methods  for  the  solution  of  stiff  ordinary 
differential  equations  leads  naturally  to  the  selection  of  strongly  stable  methods. 

The  theorem  on  convergence  (Theorem  3.  7),  which  follows  for  such  strongly 
stable  and  consistent  NLMS  integrators,  is  also  presented  in  this  section.  Three 
lemmas  needed  to  prove  the  convergence  theorem  are  developed.  Some  of  these 
proofs  are  a  modification  and  extension  of  the  proofs  used  by  Henrici  [17]  for  LMS 
methods.  The  theorem  of  A-stability  is  also  presented  in  this  section. 

3. 1.  Nonlinear  Multistep  (NLMS)  Algorithm 

For  convenience,  assume  A  to  be  a  nonsingular,  constant  matrix.  Equa¬ 
tion  (1. 1)'  can  be  written  as 

^  (e"Aty)  =  e"Atg(t,y);  y(0)  =  yQ.  (3. 1) 

Then,  integration  of  the  above  equation  over  the  interval  [t^ ,  gives 

y(t +1)  -  eiAh  y(y  +  {  n+1  eA<tn+H')  g(t\  y)  df .  (3. 2) 

t 

n 
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If  e-Atg(t,y)  is  a  slowly  varying  function,  a  simple  change  of  variable  allows 
successful  integration  by  conventional  methods.  For  g(t,y)  slowly  varying  and 
Re{X(A)}<0,  p( A)  »  1,  conventional  methods  require  a  prohibitively  small 
step  size, which  we  overcome  by  NLMS  methods.  Our  method  of  attack  is  to  express 
g(t.y)  as  a  low-order  polynomial  in  t  by  retaining  the  first  few  terms  of  the  Taylor 
series  of  g(t,y)  expanded  about  t^.  The  complete  derivation  of  this  idea  is 
properly  described  in  section  3.3,  Consistency.  When  g(t,y)  =0,  equation  (1. 1)' 
becomes  homogeneous,  i.  e. , 

y'=Ay;  y(°)  “  yQ  •  (3.3) 


Without  loss  of  generality,  we  consider  a  =  tjj  =  0.  The  solution  of  (3. 3)  is 
y(t)  =  eA*  y(0).  Consequently, 

iAh 

y'W = e  y<*n»  • 


which  is  the  rigorous  solution  in  the  absence  of  round-off  errors,  where  i  =  0, 1,2, 
. . .  ,  an  integer  index,  and  h  =  t^^  -  tjj . 

The  linear  multi-K-step  methods  take  the  general  form 


K 

Lot,  y  =  h 
i  7n+i 
i=u 


K 


i  i 


i=0 


n+i 


where  f  0  and  |  a  Q  |  +  1 0O  \  >  0  . 


o 

II 

<a 

explicit 

If  j 

,  the  method  is  j 

implicit 

The  generalization  of  (3. 4)  leads  to 


a  eAh<K-‘> 


yn+1 


K 


h  E  ♦Kl(Ah)  ®n+i  • 


i=0 


(3.4) 


(3.5) 


where  aK^0  and  |«0|  +  | M4>Ko(Ah))|  >  0. 


If 


4>KK(Ah)  0 

♦KK<Ah)  ^  ® 


the  method  is 


explicit 

1 

implicit 


Without  loss  of  generality,  we  assume  =  1  for  computational  convenience.  The 
coefficients  of  yn+j,  9n+i  depend  upon  (Ah).  We  note  that  we  need  K  starting 
values  to  proceed. 


3. 1. 1.  Starting  Procedure 

A  convergent  K-step  method  will  produce  a  uniquely  determined  sequence 

y  ,  y,  ,  . . .  ,  y  for  an  arbitrary  set  of  starting  vectors  S  ,  S.  ,  . . .  ,  5^  ,  •  In 
*o  'l  n  o  i  K-i 

practice,  we  obtain  the  starting  procedure  by  setting  the  starting  vectors  equal  to 
the  given  initial  vector  and  calculating  the  subsequent  (K-l)  vectors.  These 
starting  vectors  are  required  to  be  bounded  in  order  to  meet  the  stability  criterion. 
We  have  already  shown  that  if  the  differential  equation  satisfies  certain  require¬ 
ments,  the  solutions  of  the  difference  equation  depend  continuously  on  the  initial 
data.  Thus,  a  unique  solution  exists  and  the  numerical  solution  approaches  the 
exact  solution. 


3. 2.  Stability 

If  a  bounded  starting  procedure  yields  a  uniformly  bounded  solution  of  the 
approximating  difference  equation  to  the  differential  equation  (1. 1)  as  h-+0,  we 
say  that  the  method  is  stable.  We  can  also  describe  a  stable  method  as  follows: 
Let  SQ#  . ^K-l  denote  K  initial  vectors  ? 


||$.(h)||  <  M,  a  constant. 
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Let 


y.  =  S.(h)  ;  i  =  0,1,2,  ...  ,  K-l  . 

Then  3  a  constant  M' ,  independent  of  h,  J 

max  ||  yn^l  <M*. 

The  stability  is  determined  by  the  root  condition,  i.  e. ,  the  modulus  of  the  roots  of 
every  characteristic  polynomial  must  not  exceed  1  and  the  roots  of  modulus  1 
must  be  simple.  (The  characteristic  polynomial  is  discussed  in  this  section. )  We 
now  proceed  to  develop  the  theory  with  regard  to  the  concept  of  stability. 

If  equation  (1. 1)'  is  homogeneous,  i.  e. ,  g(t,y)  -0,  we  expect  that  the 
values  ycy  can  be  found  exactly  in  the  absence  of  round-off  error.  If  (3. 5)  is  to 
hold  when 

y<WelAV.>*  <3-6> 

then  substituting  (3. 6)  into  (3. 5)  gives 


V-  Ah(K-l)  lAh  „  , 

L»ie  e  y<y 

<=n 


(3.7) 


K 

For  y(tu)  and  (3. 7)  to  be  zero,  we  discover  that  ^  arj  must  be  0.  This  is  idea¬ 
ls 

deal  to  the  necessazy  condition  for  LMS  methods  to  be  consistent.  The  char¬ 
acteristic  polynomial,  following  Henrici  [16]  and  Dahlquist  [10],  for  LMS  methods  is 

expressed  by  K  . 

p(f)  =  ^  ai^  •  (3,s) 

i=0 


Polynomial  p(f)  is  said  to  satisfy  the  "root  condition"  if  K  roots  satisfy 
|  f  j )  <1  and  if  the  roots  satisfying  |  |  =1  have  multiplicity  1 .  Later  in  this 

section,  we  show  that  the  characteristic  polynomial  of  NLMS  methods  generalizes 


i 
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the  characteristic  polynomial  of  LMS  methods  and  that  our  particular  choice  of 
NLMS  methods  obeys  the  root  condition.  We  will  explore  the  necessary  condition 
of  stability  later  in  this  section  and  the  sufficient  condition  of  stability  in  section 
3.6.3.  Since  1  is  a  simple  root  of  the  characteristic  polynomial,  we  have 

K 

P<1)=£  *,=°’  M 

i=0 


which  is  imposed  as  one  of  the  conditions  of  consistency  for  LMS  methods. 
Formula  (3. 7)  can  be  written  as 

K 


K 

L  v 

i=0 


Ah(K-i)  Ah  i  AhK 
'  (e  f)  =e 


i=0 


i  AhK  ... 
e  P(f). 


If  all  K  roots  satisfy  the  root  condition,  p(f)  must  be  equal  to  0.  Then 

K 

£ 

1=0 


Ah(K-i)  ,  Ah  A  _ 

2^  «,  e  '  ' (e  r)  =0. 


Since  a  matrix  annihilates  its  characteristic  polynomial,  its  eigenvalues  must 
also  annihilate  the  same  polynomial.  Thus  the  above  formula  can  be  written  as 

K  aX(A)h(K-i)(eX(A)lif)i=  K 
i=0  1  i=0 

Equation  (3. 10)  is  a  set  of  n  equations  for  the  components  of  e^K  f  .  Each  equa- 

K 

tion  is  a  characteristic  polynomial  of  the  form  ^  *  =  0.  where  t  ,  which 


,<x.n= E  V(A,h^V(A)V= £  (3.10) 


XjhK 


i=0 


stands  for  ^  =  e  ^ ,  is  the  component  of  the  j-th  characteristic  polynomial 

311(1  ^i*  ,>•••»(.  ir  are  K  roots  of  the  j-th  characteristic  polynomial 

p(X,f)  with  X  =  Xj(A).  Thus,  we  have 

XhK  XhK  XT'  XhK  A 
p(X,l)  =  2J  ai<®  1  =  6  12  al=e  p(1)  =  0* 

i=0  1  i=0  1 


(3.11) 
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This  is  a  condition  of  consistency  for  LMS  methods.  Note  that  (3. 10)  generalizes 
(3.  8)  and  (3. 11)  generalizes  (3.  9).  Equation  (3. 10)  is  defined  as  the  characteristic 
polynomial  for  NLMS  methods, and  (3. 11)  is  a  necessary  condition  of  consistency 
for  NLMS  methods. 

Let  us  use  test  problem  3  (from  section  5,  Numerical  Comparisons)  as  an 
example  to  show  that  the  root  condition  is  a  necessary  condition  for  stability.  Con¬ 
sider  the  problem 


Using  the  nonlinear  multi-2-step  method  with  «0  =  4,  =  -5,  and  a 2  =  1,  we 

obtain 


with  the  initial  values 


.  2Ah  _  Ah  A  - 

46  y„-5e  y„+i+y„t2=0' 


y„  ■  a. 


y1  =  (.5791054,  -.  60958467E-02) 

and  the  step  size,  h ,  =  .  625.  Our  NLMS  characteristic  polynomial  is 

,  4  .  2Xh  _  2Xh  2Xh  2 

p(X,f)=4e  -  5  e  f  +  e  f 

=  0=e2Xh(f-l)  (f-4), 


which  has  two  simple  roots,  1  and  4.  Obviously,  this  violates  the  root  condition.  As 
we  proceed  to  solve  this  problem  numerically,  we  can  see,  from  the  computer  results 
shown  below,  that  the  divergence  becomes  evident  after  7  steps.  The  first  column  of 
the  results  is  the  iteration  index,  and  the  second  column  is  the  current  t  values.  The 
last  two  columns  show  the  calculated  numerical  values  of  two  arguments,  indicated 
by  A.  The  two  values  below  the  arguments  are  exact  solution  values,  indicated  by  T. 


\ 


J 
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1 

. 12500000+01 

,16589938+00 

.16591649+00 

-.17463091-02 

-.17464893-02 

2 

.18790000+01 

.4750o46l-01 

,47535869-01 

-.50006794-03 

-.50037757-03 

3 

.25000000+01 

.13581334-01 

,13619255-01 

-.14296139-03 

-.14336057-03 

4 

.51290000+01 

.38569114-02 

.39019817-02 

-.40599063-04 

-.41073492-04 

5 

.57500000+01 

.10657062-02 

.11179365-02 

-.11217958-04 

-.11767753-04 

6 

,43730000+01 

,26023676-03 

.52029417-03 

-.27393339-05 

-.33715176-05 

7 

.50090000+01 

.22870184-04 

.91765814-04 

-.24073868-06 

-.96595595-06 

B 

,66230000+01 

-.52686605-04 

,26291346-04 

•  55t*6lB94-06 
-.27675101-06 

9 

,62500000+01 

-.82990288-04 

,75325968-05 

.67358188-06 

-.79290492-07 

10 

.OB780000+01 

-.10158864-03 

,21581251-05 

.10693539-05 

-.22717106-07 

11 

.75000000+01 

-.11826226-03 

.01831319-06 

.12450763-05 
-, ©5085600-08 

12 

.61210000+01 

-.13609010-03 

,17714970-06 

.14325272-05 

-.18647337-08 

13 

,87500000+01 

-.15611957-03 

,50754238-07 

,16433637-05 

-.53425514-09 

14 

.93750000+01 

-.17896583-03 

,14541333-07 

,18838506-05 

-.15306666-09 

15 

,10000000+02 

-.20511784-03 

,41661615-08 

.21591349-05 

-.43854332-10 

16 

,10685000+02 

-.23508066-03 

,11936253-08 

,*4745329-05 

-.12564477-10 

17 

.11250000+02 

-.26941724-03 

,54197936-09 

.*8359707-05 

-.35997828-11 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 

A 

T 
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The  nonlinear  multi-2-step  method  used  to  solve  problem  3  is  stable  be¬ 
cause  we  chose  aQ  =  0,  =  -1,  and  Og  =  1*  so  that  the  characteristic  polynomial 

has  two  distinct  roots,  0  and  1.  Hence,  the  root  condition  is  satisfied,  and  the 
method  is  stable.  Thus,  we  have  shown  that  the  root  condition  is  a  necessary  con¬ 
dition  for  stability.  We  defer  the  proof  of  the  sufficient  condition  to  lemma  3. 6. 3. 
Thus,  the  stability  for  NLMS  methods  is  a  direct  generalization  of  the  stability  of 
I  MS  methods  (Henrici  [16]). 


3. 2. 1.  Strong  Stability 

A  particularly  advantageous  feature  of  NLMS  methods  is  that  the  strong 
stability  condition  results  when  Re  |  \(A)  |  <  0.  Since  stiff  differential  equations 
frequently  occur  when  Re  j  X  (A)  l  <  0  with  p  (A)  »  1 ,  it  is  most  important  that 
the  parasitic  growth  of  the  extraneous  solution  of  the  difference  equations  be  damped 
out.  To  ensure  this,  the  NLMS  methods  are  selected  to  be  strongly  stable.  This 
is  not  a  restriction  of  NLMS  methods  since  the  methods  are  also  applicable  when 
Re{x(A)}>0  ;  in  this  case,  the  error  growth  is  appraised  by  the  estimate  fur¬ 
nished  by  lemma  3.  6. 3. 

A  measure  of  the  growth  of  LMS  methods  is  provided  by  examining  the  solu¬ 
tion  of  the  homogeneous, constant  coefficient  difference  equation.  Consider  the 
homogeneous  equation  of  LMS  methods, 


K 


i=0 


i  *n+i 


0. 


(3. 12) 


We  recall  that  strongly  stable  solutions  of  the  LMS  methods  occur  when 

Mr)  =  £  V* » 

i=0 


I 


i 
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p(t) 

pl(f)  =7^1  has  r00tS  h’  f3’ •**’  ’ 


where 


|f2| *  |f3| » •••  »  |fK|  <  1- 


The  corresponding  homogeneous  equation  of  NLMS  methods  is 


K 


V  a,  eAh(K-‘>  y  =0. 

«  1  aH 


(3.13) 


Since 


let 


eA»>(K-l) 

»n4-i  r 


rn+i  n+i 


Then  (3. 13)  becomes 


K 


§“‘w-‘=0- 


(3.14) 


Equation  (3. 14),  of  course,  has  precisely  the  same  constant  coefficients  as  the 
linear  difference  equations  of  order  K.  In  fact, 

w  ..-eAh<,W>y 


n+i 


n+i 

Ah(K-i)  .Ah  i 
=  e  (e 

AhK  i. 

-e  fjt  3-1.2 . K. 


Strongly  stable  solutions  will  result  for 

X(A)hK 


e 


<  1. 


Note  that  these  extraneous  solutions  damp  out  extremely  fast  since  He  {  \  (A) }  <  0. 
Therefore,  we  see  that  strong  stability  implies  that  ||e^^  ^  yn+J|  *S  un^orm^ 
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bounded.  For  A  =  G.  ||  e^h^K  ^  y 


n+i 


reduces  to  ||yn+i||»  which  means  that 


||  y  +.  ||  is  uniformly  bounded;  this  is  the  stability  definition  for  LMS  methods. 


3.2.2.  A-Stability 

Another  advantageous  feature  of  NLMS  methods  is  that  when  they  are  used 
to  solve  stiff  equations, they  are  A-stable.  Dahlquist  [10]  defined  a  method  to  be 
A-stable  if  the  numerical  solution  ||  y^  ||  -*  0  asymptotically  as  n  -*  »  for 
the  differential  equation  y'=Ay  where  Re  |  \(A)}  <  0 .  If  A-stable  methods  of 
order  higher  than  2  exist,  they  do  not  belong  to  the  linear  multistep  family  since 
it  has  been  proved  by  Dahlquist  [  10]  that  an  LMS  method  of  order  higher  than  2 
cannot  be  A-stable.  However,  since  A-stability  is  a  desirable  property  when  solv¬ 
ing  stiff  equations,  it  is  preferable  that  it  be  retained  in  NLMS  methods.  We  will 
introduce  a  theorem  which  shows  that  NLMS  methods  accommodate  the  A-stability 
in  the  sense  of  Dahlquist. 

Matrix  Exponential 

We  begin  by  discussing  the  computation  of  a  matrix  exponential,  e^.  If  A 
is  a  scalar,  there  is  no  difficulty  in  computing  e^.  If  A  is  a  stable  matrix 
(Young  [33]),  i.  e. ,  Re  |X(A)[  <  0,  then  the  rational  Pade  approximation  is  also 
stable  (Varga  [31]  and  Lawson  [24]),  as  shown  by  the  following  lemma. 


Lemma  3.  2.  2.  (Pade  Lemma) 

Ah 

Denote  the  Pade  approximation  to  e  by  Pade  (Ah),  Then  for 
Re  |  x(A)  ]  <  0,  the  Pade  (Ah)  is  stable,  i.  e. ,  p  (Pade  (Ah))  <  1. 
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The  Pade  approximation  to  a  matrix  exponential,  eA,  has  many  different 
expressions,  which  can  be  found  in  references  [5]  and  [31],  In  our  present  test 
computations  following  Blue  [5],  we  use 

for  Re  |x(A)|  <  0.  Where  x  is  a  real  scalar  <0,  e*  can  be  calculated  directly 

from  an  accurate  exponential  routine.  The  Pade  approximation  is  applied  with  the 

Ah 

requirement  that  p(Ah)  <  1  when  Re  j  X  (A)  |  <  0.  The  Pade  approximation  to  e 
using  a  polynomial  of  degree  n  in  the  numerator  and  m  in  the  denominator  has  an 
error  0(hn+m+l)  as  h-*0  (Varga  [31]).  If  p(Ah)  is  not  <1,  the  accuracy  can  be 
ensured  by  the  identity  eA*1  =  (e2  mAhj2m 

Theorem  3.  2.  2.  (A-Stability  Theorem) 

When  used  to  solve  stiff  equations,  NLMS  methods  accommodate  the 
A-stability  iu  the  sense  of  Dahlquist. 

Proof:  In  the  Dahlquist  sense,  when  applying  NLMS  methods  to  the  problem 

y’  =Ay,  which  implies  g(t,y)  =0.  NLMS  methods  produce  the  approximate 

Atn  D  Ah 

solution  to  the  problem:  y^  =  e  n  yQ  =  e  yQ.  Since  the  Pade  (nAh)  is 
stable,  the  ^lim  ||yn||-*o»  ^us  establishing  the  A-stability. 

We  note  in  passing  that  the  solution  to  this  problem  is  the  principal  reason 
for  NLMS  methods  to  be  of  interest.  The  NLMS  methods  solve  this  pioblem  rigor¬ 
ously  for  every  constant  matrix  A  in  the  absence  of  round-off  error. 


V 


v 
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3. 3  Consistency 

This  section  deals  with  the  development  of  the  theory  of  consistency  in  re¬ 
lation  to  NLMS  methods.  Later  in  this  section,  we  show  that  NLMS  methods  are 
consistent  and  demonstrate  that  our  consistency  is  a  generalization  of  the  consist¬ 
ency  of  LMS  methods.  An  immediate  need  is  to  define  what  we  mean  by  consistency. 
A  method  of  (3. 5)  is  said  to  be  consistent  if 


max 

n 


K 


E", 

i=0 


Ah(K-i) 


■ h  E  *ki  <Ah)  ■ 

i=0 


n+i 


is  small  as  h  -*  0.  We  shall  show  that  the  contistency  we  will  develop  for  NLMS 
methods  actually  satisfies  our  definition. 

Problem  (1. 1)'  can  be  written  as  (3. 1),  which  is 


d  —At  —At  .  /m _ 

—  <e  y)  =  e  g(t,y);  y(0)-yo- 

Integration  of  the  above  system  over  the  interval  f  t  ,  t  .1  gives 

1  n  n+i 


y(t  .,)  -  e1Ahy(t)  ♦  /  0+1  eA^i-‘’)  ,<t'  .y)  df 


n 


-At 

which  is  our  equation  (3. 2).  If  e  g(t,y)  is  a  slowly  varying  function,  a  simple 
change  of  variable,  i.  e. , 

*  =  eA*y, 

allows  successful  integration  by  conventional  methods  such  as  the  Runge-Kutta  or 
LMS  methods.  Lawson  [24]  used  this  idea.  For  g(t,y)  slowly  varying. 
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Re  |x(A)f  <  0,  and  p(A)  »  1,  conventional  methods  require  a  prohibitively 
small  step  size,  which  we  overcome  by  NLMS  methods.  The  method  of  attack  is 
to  express  g(t,y)  as  a  low -order  polynomial  in  t,  e.g. ,  by  retaining  the  first 
few  terms  of  the  Taylor  series  of  fl(t,y)  expanded  about  t  .  For  the  moment,  let 
us  introduce  the  NLMS  operator, flf^[y(t);h] .  (The  construction  of  such  an  oper¬ 
ator  is  described  in  section  3.4,  Nonlinear  Operator.) 

Write 

<£n [y(t);h]  =  £  aieAh(K"i)  y(t+ih)  -  h  £  «>Ki(Ah)g(t  +  ih,y ) ,  (3. 15) 


i=0 


i=0 


and  introduce  the  true  operator, 

cCfyw] =  -  Ay  -  s^y) =0  • 


(3.15)’ 


For  g(t,y)CCP+1  ,  we  evaluate  the  local  discretization  error  as  follows: 


K  K 

=  £  «,  eAh(K_i)  y  (t  +  ih)  -  h  V  <t>  (Ah)  g  (t  +  ih,y) 
i=0 


K 

r 

i^ 


•  {d?  "  Ay‘  * 


(3. 16) 


The  terms  inside  {  }  of  (3. 16)  vanish  because  of  (3. 15)'.  If  we  expand  g  (t,y)  in 
a  Taylor  series  expansion  around  t  ,  we  get 

^  i0)  *<*„»  j 

9(‘.y)  -  £ - ji - -  V 

j=0  J  ’ 

&  90)  <‘n.  y(tn»  j  9<P+1>  <V  y<‘n» 

U  j!  K  n*  (p  +  1)!  '  V 


+  (t-t ) 


(t-t ) 
n 


n 

p+1 


P+1  n  (P+1)! 


(3.17) 
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In  the  above, 

Gd+1  (t  "  tn)  =  may  **  S(P  ^  y(fc))  “  1J  <*•  y (O) II 

^  t-t  CKh 

n 

is  the  modulus  of  continuity  when  g(t,y)  satisfies  the  Lipschitz  condition  only, 
0  <  d  <  1,  and  p  +  1  =  0. 

By  substituting  (3. 17)  into  (3.2).,  we  get 


.0) 


„  ,  ‘Ah  »,  f  f‘o+1  -A(*’-‘„+i)8w<t„.y(t„)) 

y(W=e  y(tn)  +  S  [  6 

J=°  *n 


n+i  n+i  ~  '2-^(t’-t)Jdt’ 

j!  n' 


.  fn+i 

+  {  e  (p+1)!  •  n 

n 


,  |U  -A(t’-W 

+  0  J  e 

t 

n 


(t*  -  t  J**1  dt’ 

Ti 

P+1 


(t'-y 


V^-y-FTirr  dt'- 


i3. 18) 


Define 


t  -A(t’-t  .) 

X  j (Ah)  =  Jnie  n  i(f-tn)'dt*. 

t 


J 


(3. 19) 


Expanding  g(t  +  ih,y)  at  t=t  and  neglecting  the  terms  modulo  (p  +  1),  we  get 


9  (t  +  ih,  y) 


®W.y>. 


j  = 


j! 


(3.20) 


By  substituting  (3. 18'  and  (3. 20)  into  (3. 16)  and  using  definition  (3. 19),  we  get 

r  m  hi  V  Ah(K-l)  lAh  „  .  .  ^  V  AhfK-i)fei  (All‘1 

T[y(t);h]  =  2^  “i e  6  yV  L  I  j!  I 

•  4  — A  4  —A  u  *  J 


i=0 


-  h 


j=0  4=0 


•(9(Pfl)(tn.y(tn»  +  *Gp+l(t-y) 


(3.21) 


) 
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Now 


,  let  us  look  at  the  last  term  of  (3. 21).  For  g(t,  y(t))  the  class  of 

polynomials  of  order  p,  the  last  term  vanishes.  For  g(t,  y(t))  £  Cr  ,  the  last 
term  will  be  O  (h^+^);  this  is  followed  by  examining  the  bound  for  the  last  term. 
The  bound  is  given  by 

P+1 


K 

E«i 

»  r\  * 


i=0 


Ah(K-i)  ±i_ 
L  CP 


(Ah) 


+  1)! 


<9 


,(P+1) 


K 


(V  ^n*)  +  SGp+l (t  ‘  n»  ||  -  E  |“il 


rP+2 


•||eA1*(K-I)l| (ii9<ptl)(‘n.y(y)ii+iiGp+1(Kh)ii)i.p+2=o(hp+2).  <3.22) 
As  h  ->  0,  this  bound  vanishes. 

For  arbitrary  g  (t  we  select  {  }  of  (3. 21)  to  be  zero,  so  that 

lim  r[y(t);h]  =  0  .  This  defines  consistency  for  NLMS  operators  and  shows  that 
h— 0 

the  NLMS  operator  is  consistent  with  the  true  operator  in  the  sense  of  Keller  [22], 

Since  the  true  operator  is  0  »  then  the  lim  T[y  (t);  h]  =  0  is  equivalent  to 

h-0 


lim  max 
h-0  n 


K 

ie 

i=0 


Ah(K-i) 


K 


yn+1  - h  E  *K,(Ah)  9 


i=0 


Ki' 


'n+i 


=  0  , 


which  is  our  definition  of  consistency.  The  selection  of  {  }  of  (3.21)  to  be  0 

enables  us  to  determine  «V,(Ah)  and  guarantees  the  consistency,  so  that  a  formal 

Ki 

analogy  with  LMS  methods  results.  Thus,  we  see  that  our  NLMS  methods  are  a 
generalization  of  LMS  methods. 
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3.4  Nonlinear  Operator 


Since  we  have  now  establish**!  the  consistency,  let  us  define 


3  _ 
i 


A  If  (Ah)  ^j+J-  r*n+i  -A(t'-twi4) 
~ =  ~  \  e 


j! 


n+1(f-t  >J  df. 
'  n' 


For  j  =  0, 1,  we  get 


S+ *  1 

] 

9‘- *!  f 


i«,  , 

e  dt'  =  e  - 1 


n+i  A(*  .  ....  .  iAh  . 

e  (t1  -  tQ)  dt’  =  (e  -  I)  -  iAh  . 


By  induction,  we  get 


(1  m+1  Q,m  (!Ah) 

**  l  ** i  "  (m  +  1)! 


then 


9 


(iAh)-fi 

Ti 


(3.  23) 


We  already  have  associated  the  nonlinear  operator  with  the  nonlinear  multi- 
step  formula  by  (3. 15).  Take  n  =  0;  then 

K  Ah(K-i)  .  K 
ae  ’y-h 

i=0  i=0 

If  we  use  formula  (3. 18)  for  and  formula  (3.  20)  for  g^  (letting  p  -*•«> )  and 
formula  (3. 23),  and  substitute  them  into  (3. 24),  we  get 


K  K 

ct  N0(t);h]  =  £  “i  eAh(K_i)  Yi  -  h  £  9^  (Ah)  .  (3. 24) 


K 


i=0 
K 


U w*i-  •»] 


«  I?  (Ah) 


^aieAh(K-l)eiAhy 
i=0  1 


(q) 


+  C^(h)g  +  C^(h)  £j T  +  •  •  •  +  Cq(h)  9  +  •  • .  • 


(3.  25) 
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Note  that  because  of  the  root  condition,  |  |  of  (3. 25)  -*0  •  Thus, 

XN[y(‘)ih]  =  E  cm9(ra). 

m=0 

where 


eAh(K-i) 


llH  (W 

1'  J  «  11 


♦Ki(Ah), 


j  =  0,1 


f  m  •  •  • 


(3.26) 


(3.27) 


We  define  the  order  p  of  NLMS  methods  to  mean  that 

CQ(h)  =  Ct(h)  =  .. .  =  Cp(h)  =0  but  Cp+1  ?  0  . 

At  this  point,  let  us  distinguish  the  orders  between  LMS  and  NLMS,  i.  e. , 

PLMS  311(1  PNLMS*  A  =  0  •  tt  is  lmPlied  that  g(t,y)  =  f  (t,y)  =  y\  Then 

(3. 25)  becomes 

<£N[y<‘>'h]  -iL[yW:h] 


CjWy 


O+i) 


a 

a=o  *; 


t)  y  +  Co(0)  y’  +  C1(0)  y"  + ...  +  Cq(0)  y(q+1)  +  . . . . 


If  we  reindex  the  coefficients  of  y^,  we  get  exactly  the  coefficients  of  the  LMS 
operator;  this  confirms  that  NLMS  methods  are  a  generalization  of  the  LMS  methods 
since  LMS  is  a  special  operator  of  NLMS.  When  A  ^  0  ,  the  two  p's  are  not  com¬ 
patible  since,  in  this  case,  the  solution  deponds  on  an  exponential  times  a  poly- 


Domlal  of  order  PNLMg .  Thee  we  have  pNLMS  =  PLMg  -  1- 
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3. 5  Formulation 
3.  5. 1  General  Formula 

By  our  selection,  we  set  {  }  of  (3. 21)  =  0  for  j  =  0, 1,  . . .  ,  p  to  give 
p-th-order  methods  ? 


f  Ak(K-i,riJi(Ah)i ,  f  <**»'*  ,4M  „ 

£v  Lin  §Tr*Ki(Ah,=o- 


(3.28) 


The  above  equation  expresses  a  K-otep,  p-th-order  method  that  consists  of  a  sys¬ 
tem  of  p  +  1  equations.  The  choice  of  K  and  p  will  determine  whether  this  system 
has  a  unique  solution,  has  many  solutions,  or  has  no  solution.  We  make  a  choice 
)  p=K,  which  ensures  the  existence  of  inverses  of  K  and  H  .  Thus,  we  can 

determine  A  ,  (Ah)  based  on  the  choice  of  a.'s,  which  can  be  written  as 
Ki  i 


i 


K 

Z  v 

i=0 


Ah(K-i) 


Aj+1l|(Ah) 


jl 


E  vAh»- 


(3.29) 


Substituting  (3.23)  into  (3.29),  we  get 
K  „Ah(K-l)  fiAh 


E  V 

i=0 


(iAh)J 

i! 


J  i=0 


1  *Ki  (Ah*  *  (3* 30) 


Because  of  the  root  condition,  (3. 30)  becomes 

f  .  PAh(K-‘>  Z  ^  -  16jP  Z  +K1  <Ah>-  * 31> 


j£=0  "•  j!  i=0 

If  we  multiply  both  sides  of  (3. 30)  by  (Ah)  and  let  ||  A  ||  -*  0,  we  get 


K,  jj+1  £  4j 

^  “i  (j  +  1)!  =  S  V  *Ki(0)  * 


(3.32) 


i=0 


Then  j  =  0  gives 


K  K 
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the  consistency  condition  for  LMS  methods  for  ▲  (Q)  =  (3„.  =  0,  •  This  confirms 

Ki  Ki 

that  NLMS  methods  are  a  generalization  of  LMS  methods. 

We  determine  ♦if.(Ah),  without  loss  of  generality ,  by  selecting  a  =  1. 

Ki  K 

In  addition,  we  require  the  condition  of  strong  stability  to  be  realized,  i.  e. , 

Re  |X(A)|  <  0.  The  (Ah)  are  determined  by  means  of  (3.31),  which  can  be 

considered  as  a  matrix  equation  for  K-step,  p-th -order  method  (K  >  1,  p  >  0).  On 

the  other  hand,  we  can  select  ♦^(Ah)  to  determine  o^'s  as  well,  but  we  choose 

not  to  do  this  since  we  would  have  to  investigate  the  strong  stability  of  the 

resulting  a.'s.  It  is  easier  to  choose  a.'s  to  be  strongly  stable  and  then 

solve  for  ^  (Ah).  An  approach  to  the  selection  of  a.'s  to  ensure  strong  stability 

Ki  i 

is  presented  below.  The  (K  x  K)  companion  matrix  of  the  characteristic  polynomial 
K 

p(f)  =  V,  ",  f  takes  the  form 
i=0  1 


K 

for  a  chosen  to  be  1.  We  require  that  p(  1)  =  0;  then  X  must  equal  0. 

K  £6  1 
Apply  Gerschgorin's  estimate  columnwise  to  obtain  the  eigenvalues  of  the  com¬ 
panion  matrix,  and  impose  the  condition  that  the  eigenvalues  lie  on  the  boundary  or 
inside  the  unit  circle.  For  the  first  (K  - 1)  columns,  we  have  the  same  estimate: 


|  X  -  0|  <  1  . 


X 


x 
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The  estimate  of  the  last  column  gives 


which  implies  that 


K-2 

I  *  -  (-"k-i’I  <-  E 

i=0 


f 


K-l 

IM  <  E  l“il 

i-0 


where  we  require  the  bound  to  be  <  1. 


Next,  we  look  for  conditions  of  that  produce  strong  stability. 


i=0 


K-l 


Consider 


where 


The  associated  (K  -  L)  x  (K  - 1)  companion  matrix  of  the  characteristic  polynomial 


I 
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Again,  applying  the  same  procedure  used  above  to  estimate  the  eigenvalues  and  to 
ensure  that  the  eigenvalues  lie  inside  the  unit  circle,  we  get 

K-3 

|X-0|  <1  and  |x  -  <-«K_2)|  <  £  l«,l  - 


which  implies  that 


Thus,  using 


K-2 

IM  <El«,l 

•  _  A  * 


i=0 


<  1  . 


=  1 


and 


K-l 

E 

i=0 

K-2 

1  K  i 

E 

E*)‘2 

j=0  1  f 

■o°‘ 

i=0 

we  can  select  atj's  satisfying  the  condition  of  strong  stability. 


3. 5. 2.  Matrix  Formula  for  $ 

K1 

In  (3.31),  i  is  a  column  index  and  j  is  a  row  index.  Let  $  be  a  vector  of 

K  - 1  or  K  elements  whose  components  are  ^(Ah);  i  =  0, 1,  . . .  ,  K  - 1  or  K. 

Then  .  (Ah)  can  be  determined  by  the  matrix  equation 
Ki 

4>  -  -K'1  H_1E  *  .  (3.33) 

Each  of  these  symbols  is  described  in  the  next  section.  For  the  system  of  equa¬ 
tions,  the  above  matrix  elements  are  all  submatrices. 
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v 


3.5.3.  Explicit  Schemes 

Equation  (3.31),  expressed  as  a  predictor,  takes  the  following  form  when 
*KK  <Ah>  = 0 : 


K 

4-i.  i 


Ah(K-i)  y'  =  _  (Ah> 


i=0 


?=0 


J2 ! 


j! 


j+1K-l 

£ 

i=0 


i3  4>Ki<Ah) 


(3.34) 


for  j  =  0, 1,  . . .  ,  p .  Now,  we  have 


(3.35) 


V 


In  the  matrix  form,  we  get 

E*  =  -HK4>. 

Thus,  $  can  be  obtained  by  (3.33).  Here, 

H  «  (Hu) 

K*es) 

e  -  (ei^) 
**  (*j,  i) 
+■  (+i.i)  • 
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(3.36) 


where 

i  =  1,  ...  i  p 

j  =  1 . K 

4  =  1,  ...  ,  K  +  l. 

3. 5. 4.  Implicit  Schemes 

If  we  use  (3. 31)  as  the  formula  for  a  NLMS  corrector,  for  (Ah)  /  O, 
we  have 


\ 


V 
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for  j  =  0, 1  . . .  ,  p.  Nov/  we  have 


(3. 37) 


This  gives  the  same  matrix  form  as  (3.36),  where 

^  =  (^mm) 

K  *  (“mm) 

1  *  (*mi) 

*  *  (*il) 

•x  *  *  (+il)  ’ 

where 

m  =  1 . p  +  1 

j2  =1,  ...  ,  K  +  1  . 


V 
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3.  G.  Lemmas 

We  now  come  to  the  most  important  theorem,  the  convergence  theorem. 
Before  this  theorem  can  be  pro\  ed,  three  lemmas  are  needed.  The  proofs  of  two 
lemmas  and  of  the  convergence  theorem  are  an  extension  and  modification  of  the 
proofs  of  LMS  methods  developed  by  Henrici  [16],  We  begin  by  establishing  the 
formula  for  discretization  error  in  the  approximate  solution  of  an  arbitrary  differ¬ 
ential  equation  by  NLMS  methods. 


3.6.1.  Lemma  3.  6.1 


Consider  the  general  characteristic  polynomial, 

K 


p(x,n  =  P(x(A),r)=E  Me  n. 

i-0 


which  satisfies  the  root  condition;  i.  e. , 


p(\,l)  =  0  V  X  =  X.(A)  . 


It  is  easily  seen  that  p(A,f)  also  equals  0  : 


0  =  ,(A,f)=X“lOAhKfi)=E“i 


Ah(K-i)  #  Ah  i  _ 
3  '(e  {•)  =0 


K 

'  .  Mats.  i 

,  ai  <e  f  >  = 

i=0  i=0 

The  following  lemma  is  a  generalization  of  Henrici  [16] 


(3.38) 


Lemma  3. 6. 1;  Define  the  scalars  7^  =  0  for  i  <  0,  where  SL  is  an  integer  index. 


Then,  3  a  set  of  bounded  scalars  {  7^  }  ? 

V  „  Ah(K-i) 

i  V(K-i) 


i=0 


I  0  =  0 

> 

0;  Jf>  0 


(3.39) 


V 
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Proof:  Write  (3  38)  as 


K 


p(A.f>  =  ,<*>-£  ",  e‘-o. 


i=0 


Ah 


where  |  =  e  f .  Consider 

#.  .  ,  ,  Ah  ..... _ ,  „ _ 

P  <£>=V  K-l6  | +  •  •  •  +a2  e  £  +aie 


Ah(K-2)  K-2  Ah(K-l)  ,K-1  AhK  .K 
>  '  '  t  j.  r,  a  '  '  £  +ae  £ 

*  o  * 


K 


-K 


Ah  -(K-l) 


'  f  (“KS  +aK-le  t '  -'5+0,  e""'"  5 


Ah(K-2)5-2  +  o  ,Ah(K-l)  e-l 
2  X 


AhK  K  ,  -lv 
+  aQ  e  )  =  £  P(t  )  . 

The  roots  of  each  row  of  p  (£)  are  the  reciprocals  of  each  corresponding  row  of 

p(£).  Note  that  e^*1  ^  0  and  p(£)  =  e^^K  p(f).  By  lemma  5.5  (Henrici  [  16]), 

the  p(f)  of  each  row  has  no  roots  outside  |  f  |  =  1  because  of  the  root  condition.  It 
#  -1 

follows  that  [p  (I)]  is  holomorphic  inside  |f|<  1  for  all  n  rows.  Using  the 

#  -1 

Maclaurin  expansion  for  each  row  of  [p  (£)]  ,  we  find 


[Apr1 


(3.40) 


i=0 


By  Cauchy's  estimate,  it  is  seen  that  all  7  are  bounded.  To  prove  (3.39),  we 

multiply  both  sides  of  (3.40)  by  p  (£)  and  equate  the  coefficients,  obtaining  (3. 39). 

Q.E.D. 


From  (3. 5),  eince  a  ^  0,  we  can  write 

K 


rn+K  nr 


K 


K-l 


-E », 


Ah(K-i) 


i=0 


W-Ev^. 

i=0 


(3.41) 


which  is  of  the  form 


y=G(y)  , 


(3.42) 
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where  y  =  yn+j^*  The  successive  iterative  form  gives 


>+l> 


=  o(y«) 


(3.43) 


for  any  initial  vector  y' 

Let  G(y)  be  defined  for  ||y||  <» ,  and  let  3  a  constant  k  9  0  <  k  <  1.  Then 
G(y)  satisfies  the  condition 

II G(y*)  -  G(y) ||  <  k  || y*  - y ||  .  (3.44) 

Using  the  definition  of  G(y),  formula  (3. 41),  and  the  fact  that  g  (t,y)  satisfies  the 
Lipschitz  condition  with  Lipschitz  constant  L,  we  see  that  condition  (3.44)  is  satis¬ 
fied  by 


h|H>KK<Ah)||  _ 

k  —  ~  L 


a 


(3.45) 


K 


.<°> 


for  sufficiently  small  h  and  for  all  ||A||  <  «°  . 

For  the  iterative  procedure  (3.43)  to  converge  for  arbitrary  initial  yv“', 
k  is  required  to  be  <1: 

h  11^  (Ah)  || 


k  <  1- 


a 


L  <  1  . 


(3.46) 


K 


Conventionally,  when  using  LMS  K-step  methods  with  ay  =  1,  we  select  h  ? 


hL*  ~  *  (<1). 


(3.47) 


Similarly,  for  NLMS  K-step  methods,  we  select  h^  to  satisfy  condition  (3.46): 


ll*KK<AhN>ll  hN  L  "  “ 


Combining  (3. 47)  and  (3. 48),  we  see  that 


V 


*KK(AVl  «kI1- 


(3.48) 


For 


WAhN> 


Q  not  too  small,  we  know  that  L*  »  L  ;  therefore,  h  »  h. 
K  N 


This  tells  us  that  we  can  choose  a  much  larger  step  size  using  NLMS  than  we  can 
choose  using  LMS. 


3.6.2.  Lemma  3. 6.  2 

Lemma  3.  6. 2:  Let  h  satisfy  the  condition  (3.46).  Then 

K 

Ell+Kl(Ah>l<“- 

i=0 

Proof:  The  4V,  (Ah)  are  linear  transformations  in  a  finite-dimensional  vector 
Ki 

space;  therefore,  they  are  completely  continuous  (Bachman  [l]).  Every  completely 
continuous  transformation  is  itself  continuous.  In  finite-dimensional  spaces,  a 
linear  transformation  is  bounded  if  and  only  if  it  is  continuous.  Hence, 

II  *K1  < Ah> II  <  ”  (3.49) 

for  i  =  0, 1, 2,  ...  ,  K ;  therefore, 

K 

EH4»Ki(Ah)|l<0°- 

i“0  r\  rv 


3.6.3.  Lemma  3.  C.  3 

The  next  lemma,  a  generalization  of  Henrici  [16],  concerns  the  error 
growth  of  NLMS  methods.  The  proof  of  this  lemma  is  preceded  by  a  list  of  neces¬ 


sary  definitions.  Let 


' V 


max|eKX^A^  h|;  0<Re{x(A)}< 
E  = 

1  ;  Re{\(A)}  <  0 

(4)  Let  z„  be  initial  guesses 
ll*J<  3Vv 


(5)  Let  \v  be  the  growth  parameters 

I X  J  <  A  V  « 


<6>  ”  =  rll<I-h*K>n-K> 

K 

(7)  X>K1<Ah>n<« 


-1 


i=0 


(8)  n  =  0,l . N. 
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Lemma  3.6.3;  The  growth  of  the  solution  (zm)  of  NLMS  satisfies  the  following 
inequality: 

\\zj  <  ff(aEQ+NA)  eDha<<>.  (3.52) 

Proof:  We  begin  by  examining  the  error  growth  of  the  NLMS  nonhomogeneous  dif¬ 
ference  tquation: 


Z 


+  « 

eAh* 

KAh 

.  +  a  e  Z 

m+K  K-l 

m+K-1 

o  m 

=  hi<* 

M 

+ 

O- 

*  .  +  . . .  $ 

\  K,m 

m+K  K-l,  m 

m+K-1  o 

+  \ 


m  * 


(3.53) 


where  m  =  n  -  K  -J? ;  i  =  0, 1 . n-K.  For  A  =  0 ,  this  reduces  to  the  LMS 

case.  For  A  ^  0 ,  we  approach  an  error  bound.  We  multiply  both  sides  of  (3. 53) 
by  7  for  .0  =  0,1,  ...  ,  n  -  K .  We  then  sum  up  each  side  and  use  formula  (3. 39)  to 


obtain: 
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LHS  -  Z  +  (a  y  __  +  ...+  a  y  )  z 

n  '  K-l  n-K  o  n-2K+l'  K-l 


+  (aK-2  6 


2Ah  ,  i  KAh 

7n-K  '  *  *  V  7n-2K+2^  *K-2  +  *  *  * 


KAh 

+a  e  y  z 
o  n-K  o 


and 


RHS  =  hK,n-KYo}  *n  +  “  {*K-l,n-K  Vo  +  *K,„-K-1  ^l}  Vl  +  '" 

+  h{*o,n-KTo  +  *l,n-K-l  T1  +  +  *K,ii-2K7k}  Xn-K  + 

“IVoVkI'.^Vk'/Vk-i  * 1 "■•••  +xoVk>- 

The  coefficients  of  t  ,  Z. . z  ,  of  the  RHS  are  functions  of  y. ,  . 

o  l  n-l  i  i,n-K-i 

If  we  apply  lemma  3.  6. 2,  the  norm  of  the  sum  of  each  term  inside  {  }  <  <j>  r. 
Equating  LHS  and  RHS  gives 

V<  >W<  >XK-2  +  -"+<  >*o-h*IC..-K»o*»  +  h<-  >Vl+- 

+  hf  >  Xo  +  Vk  +  Xn-K-1  yl  *  ■ '  ■  +  \  VK>  • 

Solving  the  above  equation  for  z  ,  we  obtain 

n 

Xn  =  <I-h,l'K.n-K>o)'1[h{  >  Xn-1  +  h  {  >  *n-2  +  - " '  +  h  {  >X<,-{<  >  XK-1 

M  >V2  +  -  +  <  >XoHVkV\,-K-iV-"+\>VK>]  '  (3-54) 


Applying  our  above  definitions  to  the  terms  of  (3. 54),  we  obtain 

K 


>XK-l  +  <  >XK-2  +  -+<  >Xol 


i(E  i-j) 

'm=0  ' 


HXn-K7o  +  Xn-K-l7l  +  •“  +  \>  7n-KH-NAr 


Er;  =  r5QE  (3.55) 


(3.56) 


n-l 


llh{  +  }*n  2  +  ...  +h{  }*oll<hr<f>  £ 


m=0 


m 


•  •  • 


(3.57) 


Then,  taking  norms  of  both  sides  of  (3. 54)  and  applying  estimates  (3. 55),  (3. 56), 
and  (3. 57),  we  find 


ll*nll  <  lla  -  h+  _K>0)“lll{r;r,E  +  rh«  ^  ||*m||  +  rN/.} 

m=0 


n-l 

m=0  7 


(3.  58) 


Let  L  =  o<t> , 

K#  =  o(QGE  +NA)  ; 
then  (3. 58)  takes  the  form 


t„iiihL#E  ii*mii  +  K#- 

m=0 


(3.  59) 


Note  that  QE<r  >  1  -*  K  >  5.  Using  mathematical  induction,  we  obtain  the  estimate 

(3  *50) 


|*m||<K#  (l  +  hL#)m 


true  for  m  =  0, 1, 
and  using 


K  - 1.  Assuming  that  (3. 60)  is  true  for  m  =  0, 1,  . , .  ,  n  -  1 


m=0 


and  the  formula 


n-l 


n 


E  ^  for  x  *  1 . 

j=0 


we  get 


XJ  <  hL1  ”£  K*<1+1,L#)“  +  K* +K#  =  K*  (1  +  hi*,"  . 
m=0  ^  ’ 


ff 

hL^ 


Therefore,  we  establish  that 


llznll<K#enhL#. 

#  # 

Substituting  the  definitions  of  K  and  L  into  (3.61),  we  obtain 


(3. 61) 


Izn||  <*(:e;  +  NA)  e 


nhc?4> 


This  is  exactly  the  inequality  (3. 52),  which  established  the  truth  for  m=n.  There¬ 
fore,  (3. 52)  holds  generally  for  m  =  0, 1, 2,  . . .  ,  N. 

Q.  E.  D. 

In  section  3. 2,  Stability,  we  mentioned  that  the  root  condition  is  a  suffi¬ 
cient  condition  for  stability.  In  this  section  we  show  that  it  is  sufficient. 

A 

We  let  f  denote  the  maximum  of  the  moduli  of  the  roots  of  the  characteris¬ 
tic  polynomial  p(f).  Let  the  first  K  initial  vectors,  ,  satisfy 

||yn||5?n3  for  n  =  0,1,2 . K-l  , 


where  ^  is  a  constant.  Let  z  be  a  set  of  starting  vectors  whose  starting  values 
satisfy 

KH-s- 

A— n 

If  we  set  z  =  f  v  and  apply  lemma  3. 6. 3  to  z  ,  we  find 
n  zn  n 

|[yn||<fn.(aE5  +  NA) 

Since  all  roots  lie  on  or  inside  the  unit  circle,  ||yn  ||  remains  bounded,  thus  estab¬ 
lishing  the  stability.  This  completes  the  proof  of  the  following  stability  theorem. 
Stability  Theorem:  A  nonlinear  multistep  method  is  strongly  stable  if  and  only  if 
its  characteristic  polynomial  satisfies  the  strong  root  condition. 

Now,  we  proceed,  by  utilizing  all  the  available  lemmas  and  definitions,  to 
prove  the  convergence  theorem. 


3.  7  Convergence  Theorem 


Theorem  3.7:  (Convergence  Theorem) 

A  strongly  stable  and  consistent  NLMS  method  is  convergent. 


X 


The  first  step  in  our  approach  is  to  estimate  the  «»-N[y<t);M  at  t  =  t^ .  In 
our  general  formula  (3.  5),  g(t,y)  is  assumed  to  belong  to  C*>+* ,  in  which  case  we 
can  directly  use  (3. 25)  to  estimate  the  nonlinear  operator.  However,  g(t,y)  may 
not  always  be  differentiable,  and,  in  this  case,  we  need  to  use  a  different  approach 
to  estimate  the  nonlinear  operator,  which  is  what  we  will  do  in  our  proof.  We  want 
to  use  the  condition  for  stability  and  consistency  to  prove  that 

lim  yn  =  yw  V  tC  [a»b]  • 

h-0 

t=t 


Note  that  the  use  of  strong  stability  gives  a  desirable  estimate  for 


t 


N[y(t);h]  since  Re  {a  (A)}  <  0.  However,  growth  estimates  may  differ,  depend¬ 
ing  on  whether  Re{x(A)}  <  0  or  Re  {A (A)}  >  0. 

Proof;  Let  y  (t)  be  the  solution  of  y’  =  Ay +g(t,y);y(to)  =  yQ 

K  K 

yn  be  the  solution  of^«e  /n+i  =  h  E  *Ki  ^n+i 

i=0  i=0 


Set 


y.  be  the  starting  values  for  j  =  0, 1,  . . .  ,  K  - 1  . 
jo 


6(h)  =  max  ||y  -  y(t  +  jh)||  , 

5  J 


and  assume  that  lim  6  (h)  =0.  We  want  to  show  that  for  any  t  C[t  ,bj  , 
h-0 
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Before  we  come  to  the  proof,  let  us  make  use  of  both  the  stability  and  the 


consistency  conditions  to  derive  an  identity  that  will  be  used  to  estimate  the  NLMS 


operator  (y  (t) ;  h] .  Formula  (3 . 3 1)  gives 


£  Ah (K-i)  <iAh>4  (Ah)j+1  ^  j  x 

e  L  ~~jt~  =  “  7'  2- 1  VAh) 

i=0  Po  i=0 


Set  j  =  0.  Simplifying  the  above  formula,  we  obtain 


This  is  a  consistency  condition  when  j  =  0. 

Define  for  t  >  0 ,  the  M.  O.  C.  (modulus  of  continuity): 

w  («)  =  max  || g  (t*  ,y)  -  g(t,y)  ||  . 
|t*-t|  <* 
t*.te[to,t] 


For  i  =  0, 1,  . . .  ,  K ,  we  can  write 


where 


g(t,y)  =g(tn,y)  +  0j  w(ih)  , 


0  4 II  <  1  and  |t  -  tn|  <  ih 


Substituting  (3.  63)  into  (3.2),  we  obtain 


(3. 62) 


(3.63) 


*  iAh.,  .  .  r*n+i  A(tn+i-t  ) 

y(W  =  e  y(tn)+J  6 


[9(tn.y)+®f  «(*)]  dt' 


Then, 

„  ,  iAh„  rVi  AW’),,,  „  .  rW /V  .. 

y(tn+1)  =  e  y<y  +  J  e  dt'  g(tn,y)+  J(  e  dt 


•  e^  w(ih)  . 


(3.64) 


If  we  apply 


>eA<W‘,>dt,_A-lll.eiAh(, 


then  (3.  64)  becomes 
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y(W =  ey (tn)  - A_1  ^  - eiAh^  g<vy>  ■ A  1  - elAh^6!  ^(ih) 

If  we  write  A  1  |l  -  eiA*H  =  0(h),  the  above  formula  becomes 

y(W =  elAhy<y  -  A_1  h  -  elAh^9(tn,y>  -  0(h)efw(ih) . 

Multiplying  both  sides  by  a  eA^^  ^  and  summing  over  i ,  we  get 


K 


Ah(K-i)  ,  . 

aie  y(W~ 


K 


Ah(K-i)  iAh 


]£  “i  e^‘*v"  *'  e*  XU 


‘  [A_1  S  “i  eAh^K_i^  (I  -  eiAh)  g(tn,y)]  -  or.  eAh^  ^  O(h)0^w(ih) .  (3 


K 

r 

i=0 


65) 


Since  the  method  is  stable,  it  must  satisfy  the  root  condition.  Therefore, 

K  K 

p(  1)  =  a  =  0  ,  which  implies  |  t  of  (3. 65)  =  e  y(t  )  X]  a-  ~  0  •  Simplify - 
i=0  1  n  i=0  1 

ing  [  ]  of  (3.  65)  and  applying  the  same  root  condition,  we  find 


i=0 


K 

z 

i=0 


9(tn.y). 


Therefore, 


f>.eAh<K-l>y(t  •1f..«Ah(K-1> 

h  *  y  n+1  u  1 


«*»•*> 


K 


Ah(K-i)  _g 


-  o(h)  £  «i  ve^w(ih) 


i=0 


And 


(3. 66) 


K  K  K 

hg  ♦Ki(Ah)fl(tn+i’y)  =  h£  ‘►Ki(Ah)fl(tn,y)  +hS  VAh)efw(ih) 


(3.67) 


i=l 


t 


I 
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Then,  (3.66)  -  (3.67)  gives 


<£n  Cv  (') = h3  -  {'a'1(E  eAh(K_1)  +AhX>Ki<Ah))  9(‘n.y)} 

-  E  (OM  “i  eAh(K_i)  +  h4»K1(Ah))  e®  o,(ih)  . 

By  the  consistency  condition  (3.  62),  j  }  -*  0  . 

ll«£[y(t);h]||  <Q#hw(ih)  , 


where 


Q#h=IXl|0(h)||  |«  |  ||eAh(K'1)||  +h||+  CAhJIl). 

i=0 

Pormula(3.6)-^N[y(t);h]  gives 
K  K 

E  eAh(K‘1)(yn+i  -  y(tn+1))  -  h  g  <t.Ki(sn+1  -  9(tn+1 .  y<tn+i»)  -Q«f  h  w(ih) , 


where 


.K 

Q  =  E  «i  eAh(K_i)  +  h*H(Ah)) 


Let  y  -  y  (t  )  =  e  .  In  view  of  the  Lipschitz  condition  for  a  ,  i.  e. , 
#n  n  n 


II 9<‘n.y<tn»  -  9(‘n.yn)ll  < L  llwy  -  yjl  • 


we  can  define 


9<t  .y(t))-9(t  ,yj  = 


g  e  for  II  •  ||  ^  0 
n  n  "  nM 

0  for  llenll=0 


so  that  we  get 


E  «Ah(K'1)  Vi  - h  E  <*b>  ViV =Qef  h-«  • 


J 
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We  now  apply  lemma  3. 6. 3  with 

e  =*  ;  j(h)  =  o 
n  n  * 

G  =la0l  +lail  +  +I°k' 

E  =  max  lies'll 
i 

g  11+KiH  Vi  iE  ll+Kill  L  =  L* 

«  =  ll(I-hc|>KK(Ah))‘1ll 

t  -  t 
n  o 
N  •— 

A  =  || O  ||  hu(Kh) 

and  obtain 

(t  -t  )<rL<t> 

||e  II  <  a(GEi(h)  +  (t  -  t  )  II O  II  0>(Kh))e  n  ° 
n  "  v  no  f 

as  h  -*  0  and  both  b  (h)  and  w  (Kh)  -*  0 . 

the  above  bound  •  -*•  0  for  every  t("[t  ,b],  establishing  the  results. 


Q.  E.  D. 


4. 


COMPUTATIONAL  CONSIDERATIONS 


In  this  section,  we  first  provide  a  general  description  of  the  different 
types  of  computational  errors  and  then  describe  how  we  treat  them  at  the  present 
test  stage.  Then  we  introduce  an  algorithm  to  compute  e^  when  A  is  a  function 
of  t. 


We  define  the  error  function  C  ,  to  be  the  first  term  of  the  initial  local 

P+1 

discretization  error  of  the  NLMS  operator.  We  perform  a  general  analysis  on 

the  error  function  Cp+1 »  which  is  dependent  on  the  characteristic  polynomial 

coefficients  and  the  characteristic  roots.  We  will  show  that,  for  LMS  methods, 

it  is  possible  to  select  values  of  a.  such  that  C  .  reaches  a  minimum.  We  will 

1  p+1 

show,  by  example,  that  these  methods  are  not  strongly  stable.  Some  interesting 
results  are  presented,  which  although  they  are  not  made  conclusive  at  this  time, 
do  provide  information  for  future  research. 


4. 1  General  Considerations 

Errors  in  computation  by  NLMS  methods  are  attributable  to  the  following 
sources: 

(1)  Input  and  output  conversion  errors 

(2)  Computational  round-off  errors 

(3)  Matrix  inversion  errors 

(4)  Pade  approximation  errors 

(5)  Local  and  global  discretization  e  rrors. 

Errors  of  types  (1)  and  (2)  depend  on  the  computing  device  and  the  soft¬ 
ware  package.  Accuracy  can  be  maintained  at  a  desired  level  by  using  double¬ 
precision  arithmetic. 
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To  minimize  the  errors,  the  matrix  inversion  package  developed  by 
Forsythe  [12]  was  used  for  all  test  problems.  If  necessary,  these  errors  can  be 
improved  further  by  using  double-  precision  arithmetic. 

The  Pade  approximation  is  stable  as  a  consequence  of  the  Pade  lemma, 
section  3.  2. 2. 

The  bounds  of  the  local  discretization  errors  can  be  estimated  by  formula 
(3.  52).  The  error  function  Cp+^  will  be  discussed  independently  in  section  4. 3. 
Since  NLMS  methods  are  strongly  stable  when  applied  to  the  solutions  of  stiff 
equations,  the  global  discretization  errors  remain  bounded  within  the  numerical 
approximation  provided  by  the  NLMS  methods.  The  growth  of  this  type  of  error 
for  all  A  can  be  estimated  by  applying  lemma  3.  6.3. 

The  round-off  errors  depend  on  the  precision  of  the  computing  device, 
which  in  turn  is  dependent  on  the  number  of  digits  used.  This  type  of  error  is 
also  dependent  on  the  number  of  operations  involved.  It  should  be  noted  that  when 
applying  NLMS  methods,  the  computation  of  eA  and  the  inversion  of  (Ah)  could 
have  been  carried  out  by  techniques  other  than  those  used  here.  To  aid  in  the 
appraisal  of  the  round-off  errors  incurred  when  obtaining  ,  we  have  pro¬ 
vided  table  1  to  show  the  number  of  operations  required  for  each  iteration  in 
terms  „f  scalar,  vector,  and  matrix  operations.  As  far  as  operational  counts 
are  concerned,  the  use  of  NLMS  methods  does  not  result  in  fewer  operations  than 
does  LMS  methods.  Albeit  the  LMS  methods  are  not  optimal  in  this  sense,  a 
submember  of  the  LMS  methods,  i.  e. ,  the  Adams  family,  does  minimize  the 
number  of  arithmetic  operations  because  it  chooses  o  =  1,  »  .  =  -1,  and  the 

remaining  o's  =0.  Of  course,  where  functions  have  operations  in  common, 


i 


t 


I 
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the  operations  should  be  calculated  beforehand  so  as  to  minimize  the  total  number 
of  calculations. 


Table  1.  Table  of  Operations 


Step  1 

Step  2 

Step  3 

Description 

Explicit 

Implicit 

Explicit 

Implicit 

Explicit 

Implicit 

I 

S 

I 

S 

I 

S 

I 

S 

I 

S 

I 

S 

(Scalar)  (Vector) 

B 

1 

1 

1 

1 

1 

(Scalar)  (Matrix) 

B 

5 

4 

7 

10 

35 

Vector  Additior 

B 

1 

2 

2 

3 

3 

4 

4 

5 

5 

6 

6 

Matrix  Addition 

B 

2 

4 

15 

22 
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(Matrix)  (Vector) 

2 

2 

3 

3 

4 

4 

5 

5 

6 

6 

7 

7 

Matrix  Multiplication 

1 

4 

5 

11 

14 

19 

Pade 

1 

* 

D 

* 

2 

* 

B 

a 

3 

a 

3 

* 

Matrix  Inversion 

1 

* 

B 

* 

1 

* 

B 

B 

1 

1 

* 

Evaluation  of  g(t,y) 

1 

1 

2 

2 

2 

2 

3 

3 

3 

1 

4 

4 

Symbols:  I  -  Initial  step 

S  -  Subsequent  steps 
*  -  Needed  if  A  =  A(t) 

4. 2  Function  A(t) 

For  Re  |  \(A(t))|  <  0  V  t,  A(t)  can  be  evaluated  by  the  Pade  approxima¬ 
tion  at  every  t^ .  Lawson  [24]  introduced  an  Algorithm  II,  a  periodic  estimation  of 


A  =  A, 

k 


af 

ay 


t=  kh 


9 


!y=y(t) 

which  can  be  adopted  in  our  evaluation  of  A(t). 

This  thesis  takes  A(t)  into  consideration,  but  in  the  test  examples  A  is 
chosen  to  be  a  real  constant  matrix.  The  evaluation  of  A(t)  will  be  included  in 
the  computation  package,  which  is  already  under  development  by  the  author. 
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4. 3  The  Error  Function  C  , 

p+1 

We  will  now  examine  the  error  function  C  ,  for  the  LMS  methods  of 

p+1 

order  p.  The  following  three  questions  arise: 

(1)  Can  we  select  a  among  the  strongly  stable  families  such  that  C  . 

i  p+1 

is  at  its  smallest  magnitude  ? 

(2)  If  the  answer  is  "yes"  to  question  (1),  is  the  Adams  family  the  optimal 
family? 

(3)  If  the  answer  is  "no"  to  question  (1),  what  is  the  choice  for  o  such 
that  Cp+1  is  at  its  smallest  magnitude  ?  Which  family  exhibits  this 
characteristic  ? 

When  answering  the  above  questions,  one  should  note  that  when  A  =  0  , 
NLMS  methods  reduce  to  LMS  methods,  if  we  multiply  both  sides  of  formula  (3. 30) 
bi  [(Ah)1*1]’  1  and  let  ||  A||  -  0  ,  we  get 

K_  J+l  .  K 


^  i  1  i 

0-n>riT  § 1  v»  • 


(4. 1) 


For  j  =  0, 1,  . . .  ,  p  -  1,  equation  (4. 1)  gives  a  formulation  of  the  LMS  methods  of 
order  p .  The  matrix  form  is  : 


/  0 


1  ...  K  \ 


\  / 


0  ±  £- 

2!  2! 


A  I 


a . 


V  s  •••  s  A-. 


/ 


K 


\  (P-1)!  (P-l)l 


h 


K0 


K1 


*  (4.2) 


V^kk/ 
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We  are  interested  in  expressing  0  (=  (0))  as  a  function  of  a. ,  and  a 

i  Ki  i  i 

as  a  function  of  the  characteristic  roots  f  .  Thus  we  need  to  perturb  the  roots  to 
determine  the  effect  on  the  error  function  C  ,  when  LMS  methods  are  used,  par- 

p+1 

ticularly  when  using  the  Adams  family  of  the  LMS  methods.  For  convenience  in 
using  formula  (4. 2),  let  us  list  the  first  few  for  LMS  methods  of  order  up 

to  3. 


Explicit  Integrator; 


For  K  =  1,  p  =  1, 


LMS  C2 
Adams  C„ 

Lt 


1 

—  a 

2  1 

1 

2  * 


For  K  -  2,  p  =  2, 


LMS  C 

o 

Adams  C 

3 


_1_ 

12 


a 


1 


1 

+  “ 
3 


a 


2 


_5_ 

12  * 


i 


I 


For  K  =  3,  p  =  3, 


(* 3,o\  (l2  ai  +  3  °2  +  4  a3  \ 


4> 


3,1 


8  4 

■"""  Q  +  ■“  Q 

12  1  3  2 


V  *3,2/  \  12  “l  +  3  “2  +  4  °3/ 

LMS 


Adams  C  =  —  . 
4  24 
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Implicit  Integrator: 
For  K  =  1,  p  =  2, 


1 

12  °1 


Adams  C3=-  . 


For  K  =  2,  p  =  3, 


A  \  /ia+i0  \ 

/  2,0\  /  12  13  2  \ 


2,1 


8  4 

—  or  +  —  o> 

12  1  3  2 


U  /  \  -  o  +  —  a  / 

^2,2/  \  12  13  2/ 


LMS  C4  =  24  °1 


Adams  C  =  -  —  . 
4  24 


"V 


For  K  =  3,  p  =  4, 


9  13 

—  a  +  —  a  +  —  a 
24  1  3  2  8  3 


19  4  9 

—  a  +  —  a  +  —  a 
24  1  3  2  8  3 


5  19 

—  —  a  +  —  a  +  —  a 

24  1  3  2  8  3 


—  a 

4  1 


3 

+  —  a 
8  3 


T  --e  r  19  1  3 

LMS  C  - - a - a - a 

5  720  1  90  2  80  3 


Adams  C  - - 

5  720 


The  resulting  C  -.In  terms  of  a.  is  summarized  in  table  2. 
p+1  i 


Let  f  be  the  roots  of  the  characteristic  polynomial  p(f)  satisfying  the 


root  condition  of  stability.  The  relationship 


n  it -t) “K=i 


(4.3) 


enables  us  to  express  as  a  function  of  .  From  (4.2)  we  can  express 

♦  (0)  as  a  function  of  a  ;  $  (0)  can  also  be  expressed  as  a  function  of  f, 

Kl  1  K1  J 


Corrector  Predictor 


according  to  formula  (4.3).  Using  f  =  1,  we  can  convert  table  2  into  table  3, 


w  1 4^ 


Using  table  3,  we  can  examine  the  change  in  C  ,  when  we  perturb  the 

p+l 


root  f.  . 
J 


Initially,  iet  us  examine  the  second-order  predictor, 


LMS  e*-s+sr2- 


The  choice  of  f  =  0  implies  that  Adams  C  =  — . 

2  O  1  2 

In  order  to  perturb  the  root  f  ,  we  write 

2 

LMSC3=^  +  E(f2  +  <>- 

It  is  easily  seen  that  when  f  +  «  =  -1 ,  C  Q  will  be  at  a  minimum.  We  vary  t  to 

perturb  f  .  When  we  keep  c  positively  small,  the  method  remains  in  the  strongly 

2 

stable  family.  As  t  —  0  and  f  -*  -1 ,  this  method  is  shifted  to  a  weakly  stable 
family.  Now  we  have  found  an  interesting  answer  to  question  (2) :  Among  LMS 
methods  of  order  p,  the  Adams  family  does  not  have  the  smallest  error  C 

p+i 

Similarly,  from  the  same  example,  we  observe  that  LMS  methods  of  order  p  that 

possess  the  smallest  error  C  1  are  not  strongly  stable.  This  can  serve  as  an 

p+l 

answer  to  question  (1).  We  now  arrive  at  the  following  conjecture:  Among  LMS 
methods  of  order  p,  the  weakly  stable  families  possess  the  smallest  error  C^, ,  • 

p+i 

The  above  study  indicates  that  there  should  exist  an  NLMS  family  that 
possesses  the  smallest  error  function  C  ,  •  This  optimal  family  is  not  yet  iden- 

p+1 

tified  and  will  not  be  identified  in  this  thesis. 


We  desire  to  make  the  best  possible  choice  of  ,  so  that  when  applying 

NLMS  methods,  a  minimum  error  function  C  ,  results.  Our  recommendation  is 

p+l 
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as  follows:  Identify  an  NLMS  family  whose  a  =  1 ,  a  =  -1 ,  a  =  a 

K  k-1  K-2  K-3 

=  . . .  =  a  =  0 .  This  family,  which  is  a  generalization  of  the  Adams  family,  we 
o 

label  GA.  We  refer  to  the  predictors  of  this  family  as  GAB  (Generalized  Adams 
Bashforth)  and  the  correctors  of  this  family  as  GAM  (Generalized  Adams  Moulton). 

We  have  made  several  numerical  investigations  with  our  test  problems  on 
different  selections  of  a., from  both  strongly  and  weakly  stable  LMS  families, 
without  noting  much  difference  in  the  results.  Without  a  thorough  round-off  error 


analysis,  one  cannot  tell  how  to  choose  a^'s  that  will  reduce  the  error  function 
C  x  •  However,  as  pointed  out  earlier,  the  Adams  fiunily  has  the  least  number 
of  operations,  and  this  will  aid  in  the  reduction  of  round-off  errors.  Until  a 


thorough  analysis  is  made  of  C  ,  .  wa  recommend  the  GA  family  as  the  represen- 

p+1 


tative  family  for  all  the  strongly  stable  NLMS  methods. 


i 


I 


5. 


NUMERICAL  COMPARISONS 


Hull  [l9]  and  Ehle  [ll]  have  thoroughly  tested  selected  existing  numerical 
methods  for  solving  initial  value  problems  in  ordinary  differential  equations;  Hull 
tested  for  both  stiff  and  nonstiff  differential  equations,  and  Ehle  tested  for  stiff 
differential  equations.  Their  work  stressed  the  requirement  that  in  order  to  com¬ 
pare  methods,  meaningful  criteria  must  be  defined.  Since  the  tests  included  in 
this  thesis  are  limited,  we  do  not  need  an  extensive  rule  for  testing  purposes. 

There  are  two  basic  reasons  why  we  do  not  intend  to  perform  extensive 

tests: 

(1)  At  this  stage,  the  principal  objective  is  to  confirm  whether  NLMS 
methods  work  effectively  on  the  selected  problems.  Some  features, 
such  as  PECm  and  double-precision  arithmetic,  are  not  yet  incorpor¬ 
ated  in  our  present  computation  package. 

(2)  Although  some  of  the  answers  that  can  be  obtained  by  using  the  Hull- 
Ehle  test  criteria  would  be  desirable,  those  answers  are  not  required 
for  our  present  purposes. 

However,  we  do  define  reasonable  te,4  criteria  that  should  result  in  a  meaningful 
comparison  of  NLMS  methods  against  the  selected  methods.  Since  each  problem  is 
different  in  nature,  we  will  define  specific  test  criteria  for  each  problem. 

To  test  NLMS  methods,  we  developed  a  program  package  for  use  with  the 
Univac  1108.  It  is  written  in  FORTRAN  V  language  and  in  single-precision  arith¬ 
metic.  Adams'  formulas  are  written  in  the  same  language  and  are  inserted  in  the 
program  wherever  needed.  Gear's  program  is  run  independently. 
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NLMS  methods  have  been  tested  extensively  on  a  number  of  stiff  problems. 
Five  stiff  problems  whose  solutions  are  known  have  been  selected  for  presentation. 
Problem  5,  selected  from  Ehle's  report,  has  four  complex  eigenvalues;  NLMS 
methods  did  very  well  on  this  problem.  Problem  6  is  presented  to  demonstrate 
how  well  an  NLMS  method  can  handle  a  nonstiff  problem;  this  problem  has  one 
eigenvalue  whose  real  part  >  0 .  The  results  of  the  comparisons  are  presented  by 
graph  or  table  or  both,  depending  on  the  need  in  each  case. 


Problem  1 


Problem  Description: 

y'  =  -lOOy  +  (1  +  t2) ;  y(0)  =  1 . 


Exact  Solution: 


2  2 
(100  t 


200t +  2) 


Problem  Parameters: 

Time  Interval:  [0,  1.  95] 

Step  Size  h  :  h  =  ;  i  =  1, 2 . 14  . 

21 


NLMS  Methods  Applied: 

Explicit  2-step. 

Compare  Against: 

Gear's  program  and  second-order  Adams-Bashforth  method. 

Comparison  Criteria: 

(1)  Relative  error  is  used  as  the  measure  of  success  when  comparing 
methods. 

(2)  Gear's  program,  Adams'  method,  and  NLMS-2-step  are  used  with 
the  same  fixed  h  =  2  over  t  C  [o,  1. 95]. 

(3)  Adams'  method  and  NLMS-2-step  are  used  with  different  h  =  —  , 

2 

i  =  1,2,  ...  ,  14  ovir  t  C  [0,  0.3],  where  the  exponential  makes  a 
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significant  contribution.  For  a  larger  step  size  >  0.  3  ,  where  the 
exponential  term  damps  out,  the  comparison  is  made  on  five  suc¬ 
cessive  computations. 

Description  of  Comparisons : 

Table  4:  Comparison  Among  Methods  of  Gear,  Adams,  and  NLMS  of 
Order  2  with  Fixed  h  =  2  14 . 

Table  5:  Comparison  Between  Methods  of  Adams  and  NLMS  of  Order  2 
for  Different  h . 

Figure  1:  t  versus  Log1Q  E  . 

Figure  2:  -Log  h  versus  Log  E . 

Ld  XU 

Eigenvalues:  -100. 

Source:  S.  Preiser  (1969). 

Remarks: 

Gear's  variable-order  technique,  as  applied  to  this  problem,  involves 

trying  different  orders  up  to  order  3.  However,  most  computations  are  car- 

-14 

riedout  with  a  second -order  stiff  method  with  an  acceptable  initial  h  =  2 

Figure  1  shows  that  for  a  fixed  small  step  size,  the  nonlinear  multi-¬ 
step  method  produces  more  accurate  results  in  terms  of  relative  error. 

Figure  2  shows  that  to  maintain  an  accuracy  of  the  order  of  10  , 

-14 

Adams'  methods  require  a  step  size  of  2  ,  whereas  the  nonlinear  multi  - 

-9 

2-step  method  can  use  a  step  size  of  2  to  maintain  the  same  accuracy. 


i 


I 


Table  4.  Comparison  Among  Methods  of  Gear  (G) 
Adams  (A),  and  NLMS  (N)  of  Order  2 
with  Fixed  h  =  2-*4 


Number  of  Steps 
t  =  nh 

Method 

Relative  Error 

G 

.2163  8476  E-03 

1 

A 

.  1055  7626  E-06 

N 

.0 

*1 

vjr 

.  7845  4426  E-02 

500 

A 

.4139  7907  E-04 

N 

.3764  9440  E-05 

G 

.3423  0498  E-02 

1,000 

A 

.  1871  0187  E-04 

N 

.  1504  4324  E-05 

G 

.  8566  0224  E-06 

5,000 

A 

.  8565  7385  E-06 

N 

.9636  4558  E-07 

G 

.6845  3327  E-06 

10,000 

A 

.  6844  9640  E-06 

N 

.4620  3507  E-06 

G 

.  1061  7635  E-05 

15, 000 

A 

.  1061  6991  E-05 

N 

. 1534  9866  E-06 

G 

.  7836  8658  E-06 

20,000 

A 

.  7836  3961  E-06 

N 

.3398  9188  E-06 

G 

.1115  4447  E-05 

25,000 

A 

.  1115  3821  E-05 

N 

.3953  2531  E-06 

G 

.  8522  7857  E-06 

30,000 

A 

.  8522  3467  E-06 

N 

.3236  3342  E-07 

Table  5.  Comparison  Between  Methods  of  Adams  (A) 
and  NLMS  (N)  of  Order  2  for  Different  h 


“  L°g2  h 

Method 

Relative  Error 

14 

A 

.4944  3871  E-06 

N 

.2317  6814  E-07 

13 

A 

.3797  0580  E-05 

N 

. 5607  4696  E-07 

12 

A 

.3049  6416  E-04 

N 

.  4306  1870  E-07 

11 

A 

.  2472  3962  E-03 

N 

.5971  7310  E-07 

10 

A 

. 2028  3254  E-02 

N 

.5312  1861  E-07 

Q 

A 

. 1685  8731  E-01 

N 

.7057  9671  E-07 

8 

A 

.  1383  0503  E-00 

N 

.  1135  2181  E-05 

A 

.  1269  6386  E+01 

7 

N 

.3016  8420  E-04 

A 

.  3223  7045  E+03 

b 

N 

.  2619  8035  E-03 

C 

A 

.  3333  8325  E+05 

0 

N 

.  1237  5420  E-02 

A 

A 

.  1516  9195  E+07 

N 

.  5426  6854  E-02 

A 

.  4206  1372  E+08 

o 

N 

.  1789  7609  E-01 

A 

.  7013  3530  E+09 

2 

N 

.3655  0522  E-01 

A 

.  7568  7259  E+10 

1 

N 

.4881  1901  E-01 

L°g  io 
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t  -  nh 


Figure  1.  t  versus  ^°g10  E 


V 
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5.  2  Problem  2 


Problem  Description: 

y*  =  -lOOy  +  (1  +  t2  -  t4);  y (0)  =  1 


Exact  Solution: 


...  99  2  .  24  \  -lOOt  /  1  2 

VlOO  “ _ 3  r _ 5  )  6  ^  ( 100  + _ 3 


24 


100  100 


100  100 


12  \  .2  4  3  1  4 

3  }  _ 2  1  100 


100 


100 


Problem  Parameters: 

Time  Interval:  [o,  20] 

Step  Size  h  :  2N(.  IE-05);  N  =  0, 1 . 

NLMS  Methods  Applied: 

Explicit  2-step. 


Compare  Against: 

Second-order  Adams-Moulton  method. 


Comparison  Criteria: 

Compare  the  results,  by  means  of  relative  errors,  after  the  first  calcu¬ 
lation  since  ihe  initial  local  discretization  and  round-off  errors  are  the 
smallest  then. 

Description  of  Comparisons: 

Table  6:  Comparison  Between  Methods  of  NLMS -2 -Step  and  Second- 
Order  Adams  Moulton. 

Figure  3 :  Log  E  versus  N . 


i 


l 
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Eigenvalues;  -100. 

Source;  S.  Preiser  (1969). 


Remarks; 

_7 

Figure  3  shows  that,  for  a  required  accuracy  of  10  ,  Adams'  step 
size  needs  to  be  chosen  ~2  (.  2  E-06);  NLMS  can  maintain  the  same  accu- 


13 

racy  using  a  step  size  ~2  (.2  E-06).  This  shows  that  1^  =  100h^  =  ||A||  h^ 


confirming  our  analysis  on  the  step  size. 


Figure  3:  Log1Q  E  versus  N 


Table  6.  Comparison  Between  Methods  of  NLMS-2-Step 
and  Second-Order  Adams  Moulton 


h  = 


(.IE -06) 
N 


Relative  Error 


NLMS  Explicit  2 


Adams  Bashforth  2 


0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


.0 

.  1490  2341  E-07 
.7451  7609  E-08 
n452  9413  E-08 
.  7455  3027  E-08 
.  1492  0056  E-07 
.  1493  8974  E-07 
.0 

.  1505  2984  E-07 
.7603  1685  E-08 
.7758  8486  E-08 
.  1615  9463  E-07 
.0 

.  1029  9568  E-07 
.7108  4830  E-07 
.7129  6095  E-06 
.  1606  4611  E-04 
.5174  4164  E-03 
.4007  2017  E-02 
.  1700  8785  E-01 
.5041  3123  E-01 


. 7450  8757  E-08 

.0 

.7451  7609  E-08 
. 7452  9413  E-08 
.7455  3027  E-08 
.0 
.0 

.2246  5323  E-07 
.5268  5443  E-07 
.4561  9011  E-06 
.3623  3823  E-05 
. 2963  6455  E-04 
. 2478  3312  E-03 
. 2167  8635  E-02 
.  2076  5414  E-01 
.2390  2048  E-00 
.3890  9526  E-K)l 
. 7090  2860  E+02 
.  2529  4871  E+03 
.  4998  2459  E-K)3 


22 


5908  2495  E-01 
1357  5798  E+01 


9109  8973  E+03 
1709  4020  E+04 
1050  7818  E+04 
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5.  3  Problem  3 

Problem  Description: 


Exact  Solution. 


Problem  Parameters: 

Time  Interval :  [0,  10] 

Step  Size  h  :  h  =  ;  i  =  0, 1,  . . .  ,  15  . 

21 

NLMS  Methods  Applied: 

Explicit  methods  of  steps  1,  2,  and  3. 


Compare  Against: 
Exact  solution. 


Comparison  Criteria; 

The  tolerance  definition  follows  Ehle;  i.  e. , 

ymaxi  =  max|||yo|| ,  max  |.y  ||  ;  k  =  0, 1,  . . .  ,  i 


Error  = 


*i+i "  *<W  V  *i> 


ymax 


i 


.-7 


Ehle's  tolerance  =  10 
Results  are  tabulated  at  t  =  10  for  different  h  by  different  NLMS 


methods  in  terms  of  above  error. 


Description  of  Comparisons 


Table  7:  Ehle  Errors  by  NLMS  Methods  for  Different  h 
(Max.  Ehle  error  =  10  8  to  10  8)  . 

Eigenvalues:  {-2,  -96 

Source:  Ehle  initial  value  problem  1  fll]  . 


Table  7.  Ehle  Errors  by  NLMS  Methods  for  Different  h 
(Max.  Ehle  error  =  10-4.  8  to  10~5.  9) 


2* 


Error 

1-Step 

2 -Step 

3 -Step 

.  1622  8589  E-ll 

.  1622  8981  E-ll 

.7368  3619  E-16 

.  1622  8589  E-ll 

.1217  2031  E-ll 

. 8113  5093  E-12 

.  1622  8981  E-ll 

.  1420  1095  E-ll 

. 1217  2423  E-ll 

.1622  9374  E-ll 

.1521  5431  E-ll 

. 1420  2273  E-ll 

.  1622  3879  E-ll 

.  1571  7888  E-li 

.1521  3075  E-ll 

.  1552  2089  E-ll 

.1597  9716  E-ll 

. 1572  6132  E-ll 

.  1327  9476  E-ll 

.  1540  0083  E-ll 

.  1597  5790  E-ll 

.  1244  0601  E-ll 

.1323  2902  E-ll 

. 1424  9808  E-ll 

.1284  1213  E-ll 

.  1241  2022  E-ll 

.1274  2658  E-ll 

.  1451  7434  E-ll 

.1283  5312  E-ll 

.  1242  0126  E-ll 

.  1637  2653  E-ll 

.  1462  6200  E-ll 

.1350  7056  E-ll 

.  1646  6864  E-ll 

.  1646  5294  E-ll 

. 1565  7965  E-ll 

.  1682  9968  E-ll 

.1682  0547  E-ll 

.  1682  0547  E-ll 

.  1757  1092  E-ll 

.  1757  1092  E-ll 

.  1756  2456  E-ll 

.  1857  5613  E-ll 

.1857  5613  E-ll 

.  1857  5613  E-ll 

Problem  4 


Problem  Description: 


Exact  Solution: 


/  -0.  It  -50t\ 
/  e  +  e  \ 


y(t)  = 


-50t 


-120t 


e 


-50t  / 
+  e  / 


Problem  Parameters: 

Time  Interval:  [0,  10] 

Step  Size  h  :  0.  01  and  0. 2. 


NLMS  Methods  Applied: 

Explicit  methods  of  steps  1,  2,  and  3. 


Compare  Against: 
Trapezoidal  rule. 


Comparison  Criteria: 

(1)  The  error  definition  follows  Seinfeld;  i.  e. , 

yn  “  y  (V 

Error  R  - — —  for  each  component. 

y<‘n> 

(2)  The  computation  time  is  also  tabulated. 
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Description  of  Comparisons: 

Table  8:  Relative  Error  Comparison  Between  NLMS  Methods  and 
Trapezoidal  Rule. 

Eigenvalues:  |-0. 1,  -50,  -120 } . 

Source:  Seinfeld,  Lapidus,  and  Hwang  [29] . 

Remarks : 

Note  that  NLMS  methods  of  different  steps  produce  errors  of  the  same 
order  of  magnitude.  This  is  expected  because  we  designed  the  methods  to 
solve  the  problem  y=  Ay  effectively. 


Table  8.  Relative  Error  Comparison  Betweu .  NLMS  Methods 

and  Trapezoidal  Rule 


Method 

h 

R 

1 

R2 

B3 

Time 

t  =  0.4 

t  =  10 

• 

o 

h 

+» 

t=  0.4 

(sec) 

Trapezoidal  Rule 

.20 

l.OxlO-3 

2.7xl0~4 

6. 5  x 107 

5 

1.3x10 

1.3 

NLMS-l-Step 

.01 

1.9  xio"5 

4,  7  x 10~4 

-6 

2. 3  x  10 

-6 

2.  5  x  10 

<1 

NLMS-2-Step 

.01 

1.9  xlO-5 

4.7  xlO"4 

2.3  xlO"6 

2.  6  xio"6 

<1 

NLMS -3 -Step 

.01 

1.  9  xio"5 

4.  7  * 10~4 

2.4  xlO-6 

2.  6  xlO”6 

<2 

NLMS-l-Step 

.20 

1.9  xio  '5 

4.  5  x 10"4 

4.  0  x 10~6 

4.0x10“ 

<1 
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5.  5  Problem  5 

Problem  Description: 


V  = 


/-104 

io3 

0 

°\ 

A 

-104 

0 

1 

0 

0 

1 

0 

-50 

10  j 
1 

l  0 

0 

-10 

-50/ 

y(0)  = 


A 

i 

l 

v/ 


Exact  Solution: 


yx(t)  =  e 
y2(t)  =  c 


-10  t 

-io4t 


cos(10  t)  +  sin(10  t) 
cos(10^  t)  -  sin(10^  t) 


-50 1 

yg(t)  =  e  [cos  (lOt)  +  sin(10t)] 

-50 1 

y4(t)  =  e  [cos (lOt)  -  sin(10t)]. 

Problem  Parameters: 

Time  Interval:  [0,  3J 

Step  Size  h  ;  3/2^  ;  JL  =  0, 1 . 15 . 

NLMS  Methods  Applied: 

Explicit  method  of  step  1. 

Compare  Against: 

Exact  solution. 

Comparison  Criteria: 

The  tolerance  definition  follows  Ehle;  see  problem  3.  Results  are  tabu¬ 
lated  at  t  =  .  91552374 E-01  for  different  h  in  terms  of  above  error. 
Description  of  Comparisons: 

Table  9:  Largest  Ehle  Errors  by  NLMS-l-Step  for  Different  h 

—2  6  —3  1 

(Max.  Ehle  error  =  10  '  to  10  ‘  )  . 
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Eigenvalues:  j-50  +  10i,  -104  ^  103  i  [  . 

Source:  Ehle  initial  value  problem  3  [ll]  . 

Remarks : 

Note  that  table  9  indicates  increased  errors  with  decreased  step  sizes; 
this  is  probably  due  to  round-off  errors.  In  addition,  it  is  not  surprising  that, 
for  the  largest  step  sizes,  errors  become  zero  since  NLMS  methods  were  de¬ 
signed  to  solve  this  type  of  problem  exactly  in  the  absence  of  round-off  errors. 


Table  9.  Largest  Ehle  Errors  by  NLMS- 
1-Step  for  Different  h 
(Max.  Ehle  error  =  10-3* 6  to  10-3, 
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5.  6  Problem  6 


Problem  Description: 


Problem  Parameters; 

Time  Interval:  [0,  0.  082] 

Step  Size  h  :  1/2®  . 

NLMS  Methods  Applied: 

Explicit  NLMS-2-step  method. 

Compare  Against: 

Exact  solution. 

Comparison  Criteria: 

Using  vector  ||*  ||  ,  compare  results  against  exact  solutions. 

CO 

Description  of  Comparisons: 

Table  10:  Table  of  Numerical  Results  and  Exact  Solutions. 
Eigenvalues:  {-10,  1}  . 


Source: 


S.  Preiser  (1969). 
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Remarks: 

NLMS  methods  were  designed  to  be  effective  for  g(t,  y)  belonging  to 
the  class  of  slowly  varying,  low-order  polynomials;  therefore,  it  is  not 
surprising  that  the  NLMS-2-step  method  produces  accurate  results. 


Table  10.  Table  of  Numerical  Results  and  Exact  Solutions 
(N:  NLMS-2-Step;  T:  Exact  Solution) 


t 

Method 

Solution 

Relative  Error 

yx(t) 

y2(t) 

.  7812  5000  E-02 

N 

T 

.  1015  6862  E+01 
.  1015  6862  E+01 

.  1015  6862  E+01 
.  1015  6862  E+01 

.  1467  1029  E-07 

.  2343  7500  E-01 

N 

T 

.  1047  4286  E+01 
.  i'4r  4286  E+01 

.  1047  4286  E+01 
.  1047  4286  E+01 

0 

.4296  8750  E-01 

N 

T 

.  1087  8105  E+01 
.  1087  8105  E+01 

.  1087  8106  E+01 
.  1087  8105  E+01 

.  1369  8305  E-07 

.  6250  0000  E-01 

N 

T 

.  1128  9889  E+01 
.  1128  9889  E+01 

.  1128  9889  E+01 
.  1128  9889  E+01 

.  2639  7356  E-07 

.  8203  1250  E-01 

N 

T 

.  1170  9794  E+01 
.  1170  9795  E+01 

. 1170  9794  E+01 

. 1170  9795  E+01 

— 

. 5090  1528  E-07 

v 
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FUTURE  RELATED  RESEARCH 


The  following  areas  of  future  research  are  desirable: 

(1)  Develop  a  package  of  NLMS  computer  programs  with  the  following 
features : 

(a)  The  inclusion  of  I  MS  methods 

(b)  Freedom  to  select  a 

(c)  Built-in  PCm  procedure  with  variable  step  size 

(d)  Double-precision  option 

(e)  Evaluation  of  A(t)  at  t  =  t^  and  application  of  Pade  approximation 
to  Aftj)  . 

With  feature  (a),  NLMS  methods  can  solve  nonstiff  equations. 

With  feature  (e),  NLMS  methods  can  solve  nonlinear  equations. 

(2)  Perform  a  thorough  round-off  error  analysis  for  NLMS  methods. 


(3)  Develop  theory  to  appraise  the  error  function 


r  so  that  for  a  cer- 

V* 


tain  choice  of  the  combination  of  characteristic  roots,  one  can  minimize 


the  error  function  C  • 

^p+1 
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7. 


CONCLUSIONS 


We  have  developed  a  family  of  strongly  and  asymptotically  stable  NLMS 
methods  for  solving  stiff  equations.  The  advantages  of  NLMS  methods  have  already 
been  demonstrated  in  various  sections  of  this  thesis.  The  following  conclusions 
can  be  drawn  from  the  woi-k  in  the  preceding  sections: 

(1)  Round-off  error  consideration  is  absent  in  tnis  thesis.  While  it  is  pos¬ 
sible  that  round-off  errors  will  be  large  for  some  problems,  all  test 
results  produced  by  NLMS  methods  using  single-precision  arithmetic 
seem  negligibly  affected  by  the  round-off  error  encountered. 

(2)  It  is  theoretically  true  that  on  a  per-step  basis,  higher  order  methods 
give  better  accuracy.  However,  it  appears  that  low-order  NLMS 
methods  (which  require  fewer  steps)  maintain  adequate  accuracy,  so 
that  we  believe  it  is  not  necessary  to  employ  high-order  NLMS  methods. 
However,  if  desired,  high-order  NLMS  methods  can  be  derived  easily. 
The  use  of  high-order  NLMS  methods  can  increase  computation  time 
and  round-off  error  and,  in  addition,  would  require  an  accurate  Pade 
approximation  and  matrix  inversion. 


(3)  The  findings  from  the  analysis  of  C 


p+l 


present  some  preliminary  in¬ 


formation  for  developing  criteria  for  selecting  the  characteristic  roots 


needed  to  minimize  Cp+^ .  Dahlquist  [10]  mentioned  that,  for  a  cer¬ 
tain  class  of  equations,  it  would  be  desirable  to  have  an  A-stable  method 


with  a  small  truncation  error. 


Once  the  minimization  of 


is  ob¬ 


tained  for  NLMS  methods,  the  NLMS  methods  seem  to  provide  this. 
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(4)  We  have  demonstrated  that  NLMS  methods  avoid  the  use  of  small  step 
size  when  solving  stiff  equations.  Not  only  are  low-order  NLMS 
methods  good  predictors  and  correctors,  they  can  also  be  used  as 
starters. 

(5)  After  he  performed  a  series  of  tests  on  a  set  of  selected  stiff  equations 
with  a  set  of  existing  stiff  methods,  Ehle  [11]  concluded  that  none  of 
the  methods  used  was,  by  itself,  suitable  for  solving  the  entire  collec¬ 
tion  of  the  selected  problems.  Means  have  not  yet  been  developed  for 
NLMS  methods  to  handle  A(t);  otherwise,  NLMS  methods  have  been 
applied  to  solve  some  of  Ehle's  problems  without  a  single  failure.  The 
NLMS  methods  consistently  produced  acceptably  accurate  results.  It  is 
felt  that  when  the  A(t)  feature  is  included,  NLMS  methods  will  be  able 
to  handle  a  much  wider  class  of  stiff  problems.  At  this  stage,  we 
can  say  that  NLMS  methods  are  effective  for  equations  whose  solu¬ 
tions  are  asymptotically  stable.  Indeed,  for  equations  of  the  form 

^  =  Ay  +  t)  (where  (P(t)  is  a  low -order  polynomial  in  t)  in 
dt  •  v  w  n 

the  absence  of  round-off  errors,  the  technique  is  exact  even  for  non- 
asymptotically  stable  ordinary  differential  equations. 


8.  APPENDIX  -UNIVAC  1108  FORTRAN  V  COMPUTER  PROGRAMS 

This  section  lists  the  computer  programs,  written  by  the  author,  that  are 
used  to  perform  numerical  experiments  by  NLMS  methods.  The  programs  are 
written  in  FORTRAN  V  language,  in  independent  subroutines  of  the  CALL  type,  for 
use  on  the  Univac  1108  computer.  Compilation  was  done  by  the  EXEC  8  system. 

At  the  present  time,  numerical  calculations  are  performed  in  single-precision 
arithmetic. 

The  variable-step-size  technique  was  not  applied  in  testing  the  selected 
stiff  problems.  The  technique  will,  however,  be  incorporated  in  the  computation 
package.  This  package  will  have  an  executive  program  to  control  the  variable 
step  size  and  will  call  the  listed  subroutines  as  required. 


77 


78 


c  ****  NONLINEAR  multi-k-step  methods 

SUBROUT INE  NLMSX ( KSTEP » HSIZE , Y , NORD , ALPHA »A»YN»T» INDEX » IS, HOLD ) 

C  ****  INl)EX=0  CALLS  FOR  PREDICTOR,  OTHERWISE,  CORRECTOR 
DIMENSION  Y (4,3b) , ALPHA (1),A(35,35),VN(1),T(1) 

DIMENSION  Uhl (4, 35,35), PHI (4, 35, 35), Pi (35, 3b) 

DIMENSION  UixlT ( 35, 3b) , AH (35, 35)  ,AH2(35»35)  ,AH3(3b,35)  ,Ah4(35,35) 
DIMENSION  c.AH(3b,35),E2AH(3b,35),E3AH(35,35) 

IF(HSIZE-HOLD)  11,10,11 

10  IF(IS  ,GT,  1)  60  TO  (100,200,300),  KSTEP 

11  CONTINUE 

DO  1  1=1, NORD 
00  ?  v'=l,NOKU 
P1(I,J)=0.0 
UNIT(I*J)=0.0 
£AH( I , J)=0*0 
E2AH(I»J)=0.0 
E3Ah(I»J)=0.0 
AH2(I*J)=0.0 
AH3( 1, J)=0«  0 
AH4 ( 1 , J) =0 • 0 

2  CONTINUE 
UNIT ( I , I ) =1 • 0 

IF  ( INDEX  .EO.  0)  AH3 ( I , I ) =1 , 0 
AH4( I , I )=1 • 0 
1  CONTINUE 
DO  3  1=1*4 
DO  4  J=1,NUKU 
DO  5  K=1 ,  NORD 
PHI ( I  * J*K)=0,0 
<2Hl  (1*J»K)=0,0 

5  CONTINUE 
4  CONTINUE 

3  CONTINUE 

DO  6  1=1, NORD 
00  7  J=l,NOHu 
AH( I , J)=HS1ZE*A ( I • J) 

7  CONTINUE 

6  CONTINUE 

60  TO  (100,200,300)  ,  KSTEP 
C  ****  MULTU-l-STEP 

100  continue 

CALL  STEP1  ( A ,  NORD ,  KSTEP ,  HSIZE ,  Y ,  UNI  T ,  PI ,  AH ,  EAH ,  PH  I ,  ALPHA ,  YN ,  INDEX , 
*T,IS»HOLO) 

RETURN 

C  ****  MUlTI-2-STEP 
200  CONTINUE 

CALL  STEP2  ( A ,  NORD , KSTEP , HSI ZE ,  V ,  UNI  T ,  PI ,  AH ,  AH2 , EAH , E2 AH , PHI ,  ALPHA , 
$YN»lNDEX,T, IS, HOLD) 

RETURN 

C  ****  MUlTI-3-STEP 
300  CONTINUE 

CALL  STEP3  ( A ,  NORD ,  KSTEP ,  HSI ZE ,  Y ,  UN  I T ,  AH ,  AH2  •  AH3 ,  AH4 ,  E  AH ,  E2  AH ,  E3AH , 
SPHI.QHI,  ALPHA,  YN,  INDEX,  T,  IS,  HOLD) 

RETURN 

END 
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NONLINEAR  MULTl-l-STEP 

SUBROUTINE  STEP1(A»N,KSTEP»H»Y,UNIT»P1,AH,EAH,PhI» ALPHA# YN* INDEX ;T 
St  IS. HOLD) 

DIMENSION  Pl(35»35),AH(35»3b>,EAH{35»35),UNlT(35»35)»Y(4,35) 
DIMENSION  A(35,35),PHI(4»35,35),G(3b),ALPHA(l)»YN(l),T(l) 

DIMENSION  Ah2(35,35) ,E2AH(35,35) 

DO  132  151#N 
00  133  J=X.N 
Pl(ItJ)=0.0 
133  CONTINUE 

IF(lS  .to.  1)  P1(I»I)=1.0 
YN(I>=0.0 
PHl(2»I,l>=0.0 
PHl(3»I»l>=0.0 
132  CONTINUE 

IF(H-HOlO)  141.140.141 

140  IFdS.GT.l  .AND.  INDEX. EQ.O)  60  TO  131 
IFdS.GT.l  .AND.  INUEX.tQ.l)  GO  TO  170 

141  CONTINUE 

C  #*♦#  pi=(ah)  inverse 

IF (N  -1)  121.120,121 

120  Pl(l,l)=l./AH(l,i) 

GO  TO  122 

121  CALL  INVERT ( AM, N, PI) 

122  CONTINUE 

C  ****  EAh=PAUE(EXP(AH)  ) 

IF(N  -1)  124.123,124 

123  EAH(1»1)=EXP(A(1.1)*H) 

GO  TO  125 

124  CALL  PADE(A.H.EAH.N) 

125  CONTINUE 

DO  103  1=1. N 
DO  104  J=1,N 

C  ****  AH(I.J)=AuPhA-0*EXP(AH)*I 

AH(It J)=ALPhA(l)*EAH(I,J)4UNlT(I»J) 

104  CONTINUE 
103  CONTINUE 

C  ****  PREDICTOR 

IF (INDEX  .NE.  0)  GO  TO  150 
C  ****  TO  COMPUTE.  PHI  (1,0) 

DO  105  I=1»N 
DO  lC6  J=1,N 
DO  107  K=1»N 

PHi (1 , I » J) =PhI (1 ♦ I » J)  "Pi  (I,K)*AH(K,D) 

107  CONTINUE 
106  CONTINUE 

105  CONTINUE 

131  CALL  GFN(G»H,N,Y,KSTEP,T»A) 
c  ****  Y(N*1)=YN 

DO  108  1=1, N 
DO  109  J=1,N 

YN( I )=YN( I ) +h*PHI ( 1 , 1 , Jt  *G( J) 

YN(1)=YN(1)-ALPHA(1)*EAH(I,J)*Y(1,J) 

109  CONTINUE 
10B  CONTINUE 


KtTuRN 

c  ****  cohRectok 

C  ****  TO  COMPUTE.  PHJ(1,0)  *  (111) 

150  CONTINUE 

IF  ( fxt  -  1)  15*:  *  151 » 152 

151  aH2 ( 1 » 1 )  -PI  il*l) *P1 (1*1) 

60  TO  153 

15*;  uo  154  I  =  1*N 
UO  155  J=1>N 
UO  156  K=1»H 

AH2(  I » J)=Artt;  1 1  #U)+Pi  (I#K)*Pl(K*  J) 

15b  CONTINUE. 

155  CONTINUE 
154  continue 
153  CONT 1NUE 
uo  160  1=1 » r. 

UO  161  J=l#h 
UO  162  K51>N 

phi  ( l » 1 1  u  >  =Hn  Kiihj)  ♦alpha  {  i  >  *h*  a  ( i »  k ) *e ah ( k .  j  j 

162  CONTINUE 

PHl(l>I.J)=HHI(l»I»  J)-Ah(l»J) 

E2AhU»J)=mh(1»U)+H*A(I,J) 

161  CONTINUE 
160  CONTINUl 
170  k2=aSTEP+1 
uO  163  i  =  l»N 
UO  164  K=l»r\^ 

CALl  GFN(6»H»N,Y»K»T#A) 

UO  lo5  J=1»N 

IFU  .EQ.  1)  PHI(3»I#l)=PHl(3fI»l)+PHI(l,I,u)*G(J) 

1F(K  .tW.  <:)  PHI (2  » 1 » 1 )  =PHl  (2»  1  *  1 )  +l2AH(  I » U)  *G  ( j) 

165  CONTINUl 

Pril  (3f  I » 1) -Phi  ( 0»  I » 1  )+PHl  l2»  I » 1 ) 

164  CONTINUE 

163  CONTINUE 

UO  166  1  =  111, 

UO  167  *=1»N 

YN ( 1 ) = YN ( 1 > -h* AN2 <i»K)*PHI(3*K»l> -AlPHA ( 1 ) *lAH (I*K)*Y(1,K) 
167  CONTINUE 

166  CONT INUl 
RETuKN 
ENU 
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NONLINEAR  MULTI-2-STEP 

SUBROUTINE  STEP2< A.N.KSTEP.H.  Y.UNIT.Pl. AH, AH2.EaH.E2AH.PHI, ALPHA, Y 
SN'InDEX'T'IS'HOlD) 

DIMENSION  PI (35,55) , AH (35. 3b)  »,iH2  (35,35)  ,EAH(35,35)  ,E2AH(35.35) 
DIMENSION  UNIT (3b>35)  >  Y(4, 3b)  >A(35>35)  >PHI  (4,35#  35) ,  AlPHA(I) 
DIMENSION  n(3>3>35>35)  >  AH3(35>35)  , 6(35)  ,YN(1) ,T (1 ) ,QHI (4*35. 35) 
DIMENSION  t3AH<35,35) 

DO  240  I=1»N 
DO  241  J=1»N 
P1(I»U)=0,0 
241  CONTINUE 
VN(I)=0,0 
PHI(4»20>I)=0.0 
IFdS  ,E«.  1)  P1(I>I)=1.0 
240  CONTINUE 

IF(IS  .6T.  1)  GO  TO  21 2 
C  ****  £AH=PAU£(fcXP(AH) ) 

C  ****  E2 AH=P aDE  ( EXP ( 2AH ) ) 

IF (N  -1)  206,207,208 

207  EAH<1»1>=EXP(AH(1,1) ) 

E2Ah  ( 1  >  1 )  =EAH  ( 1  >  1 )  AH  (1*1) 

GO  TO  212 

208  CALL  PADE(A>H>EAHiN) 

H2=H+H 

CALL  PAOE  ( A  >  h2  >  1 2 AH  >  N ) 

212  CONTINUE 

IF  (INDEX  .GT.  0)  GO  TO  260 
C  PREDICTOR 

IF(H-HOLU)  251>250>251 

250  IF(IS  .GT.  1)  GO  TO  234 

251  CONTINUE 

c  MM  AH2(I>J)  =  (AH)  INVERSE  SQUARE 
IF(N  -  1)  202>201>202 

201  AH2(1,1)=1./AH(1,1)**2 
GO  TO  206 

202  CALL  INVERT ( AH.N.P1) 

DO  203  1=1, N 

DO  204  J=1»N 
DO  205  K=1>N 

AH2(I»J)=AH2(I,  J)+Pl(InC)*Pl(K>J) 

205  CONTINUE 
204  CONUNUt 

203  CONTINUE 

206  CONTINUE 

C  ***«  EAh=PAUE(c.XP(AH)) 

C  ****  COMPUTE  PHI (2> 0) >  (2.1) 

DO  213  1=1. N 
DO  214  J=1.N 

PI  ( I  >  J )  = ALPHA  ( 1 )  *E2AH  ( I .  J )  ♦  ALPHA  ( 2 )  *E  AH  ( I .  J )  +UNI T  ( I » J ) 

214  CONTINUE 

213  CONTINUE 

DO  215  1=1. N 
DO  216  J=1,N 
00  217  K=1»N 

PHI(1’I.J)=PhIU.I.J)+ALPHA(1)*AH(I.K)*E2AH  K.J) 
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PHl(2*I#J)=PHl(2»I»J)«>ALPHA<2)*AH(I,K)*EAH(K»d) 

217  CONTINUE 
21b  CONTINUE 
215  CONTINUE 

DO  210  I=1»N 
00  219  J=1*N 

PHl<l*I.U>=PHl(i*I»J)-Pl(I* J)-AH(I»J) 

PHI  {2* I *0)=PHl  (2»  I » J)  4-Pl  { I’J)+2«*AH(IfJ) 

219  CONTINUE 
Pl(l#I)=Q.O 

210  CONTINUE 
234  CONTINUE 

00  220  K=1*KSTEP 
CALL  &FN(G*H*N»Y*K»T*A) 

00  221  I=1*N 
00  222  J=1*N 

YN(I)=YN<1)+PHI IK. I*J)*G(J)*H 
222  CONTINUE 
221  CONTINUE 

220  CONTINUE 

00  230  I=1*N 
DO  231  J=1*N 

Pl(l*I)=Pl(l*I)4>AH2(I*  J)*YN(J) 

231  CONTINUE 
230  CONTINUE 

00  232  I=1*N 
00  233  J=1*N 

PHK4»20»I)=PHI(4»20»I)-ALPHA(2)*EAH(I»J)*Y(2»J)-ALPHA(1)*E2AH(I,J 

S)*Y(1*J) 

233  CONTINUE 

YN ( I ) =P J ( 1  *  I ) +PHI (4*20*1) 

232  CONTINUE 
RETURN 

C  ****  CORRECTOR 

260  CONTINUE 

IF(IS  ,GT.  1)  GO  TO  204 
IF(N  -1)  2o2»261«262 

261  AH2  (1*1) =  An  (1*1)  *  AH  (1*1) 

AH3(  1 *l)=lt/( AH2 (1*1) *AH (1*1)) 

GO  TO  2b3 

262  00  264  1=1 «N 
00  265  J=1*N 
00  266  K=1*N 

AH2(  I  * J)=AH2 ( I  * J) ♦AH( I *K)*AH(Ki J) 

266  CONTINUE 
265  CONTINUE 
264  CONTINUE 

00  267  1=1 *N 
00  260  J=1*N 
00  269  K=1*N 

E3AH  ( I  • J ) =E3AH ( 1  *  J ) ♦ AH2 ( I  *  K ) # AH ( K  *  0 ) 

269  CONTINUE 
260  CONTINUE 

267  CONTINUE 

CALL  INVERT (E3AH*N*AH3) 

263  CONTINUE 

204  IF (IS  .GT.  1)  GO  TO  279 
C  ****  COMPUTE  PHI (2*0) *  (2*1)  S  (2*2) 
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DO  270  I=1#N 
DO  271  J=1#N 

W(ltl»X»J) SUN IT(I»J)-ltb*AH(I»J) + AH2 ( I  *  J ) 
W(l,2»I,J)=UNIT(I,J)-0#b*AH(I,J) 

W(l*3»I,U)=UNlT(I»J)+0.b*AH(I,J) 

N  (2»  1 » I » J)=-2#* (UNIT ( I » J)«Ah( I » J) ) 

W (2#2* I  * J) =-2 • *UNIT ( I  * J) *AH2 ( I » J) 

W (2#  3»  I »U) =-2« * (UNIT ( I » J) +Ah ( I »U) ) 

W(3,l»I.U)=Ml>2»I.U) 

W(3,2»I.J)=M1»3#I*J) 

W  ( 3 » 3 » I » J ) =UN I T ( 1 1 J i +1 • b* AH ( I • J ) ♦ AH2  ( X » J ) 

271  CONTINUE 
270  CONTINUE 

DO  263  K=l»3 
DO  272  I=1»N 
DO  273  J=1»N 
00  274  L=1’N 

Gi  I  iK» I » J)=UhI  (K#  I »U)+ALPHA ( 1  )*W (Kr  1'  I »L)*E2AH(L’ J)+ALPHA (2)*W(Ki2 

*,lfL)*£AHU,J) 

274  CONTINUE 

QHI(K»I,J)50HX(K.I#  J)tW(K*3, 1,J) 

273  CONTINUE 

272  CONTINUE 
263  CONTINUE 

DO  275  L=lt  3 
DO  276  1=1»N 
00  277  JS1«N 
DO  276  K=1*N 

PHl(Lf  I#J)=PhI<L»I#J)-AH3<I»K>*QHI(L*K»J> 

278  CONTINUE 
277  CONTlNUt 
276  CONTINUE 
27b  CONTINUE 

279  CONTINUE 

DO  260  K=l»3 
CALL  GFN(Gih,N,T,K»T»AJ 
DO  261  1=1 »N 
DO  262  J=1#N 

YN(I)=YN(I)+H*PMI(K»I»J)*6(J) 

IF(K  .EG*  3)  YN(I)=YN(I)-ALPHA(1)*E2AH(I# J)*Y(1# J)-ALPHA(2)*EAH(I# 
SJ)*V(2»J) 

262  CONTINUE 
261  CONTINUE 

280  CONTINUE 
END 
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NONLINEAR  MULT I -3-STEP 

SUBROUTINE  STEP3(A'N'KSTEP'H*Y'UNIT'AH'AH2'AH3'AH4'iiAH'E2AH'E3AH»P 
SHI »QHI » ALPHA »  yn*  index » t»  is» holo ) 

DIMENSION  AH(35»35)#AH2(35»35)»AH3(35,35)*AH4(35*35) 

DIMENSION  UNIT  (35 *35)  «£AH(35»35)  *E2AH( 35*35)  *E3aH(35*35) 

DIMENSION  A (35, 35)  ,Y(4;>35)  ,PHl  (4,35,35)  .QHI  (4,35,35) 

DIMENSION  w(4,4,35»35)  ,6(35)  ,ALPHA(1)  ,YN(1> » T(I) 

KUP=KST£P 

IP (INDEX  ,EQ.  1)  KUPsKSTEP+i 
DO  320  1=1, N 
YN(I)=0.0 
320  CONTINUE 

IF(H-HOLO)  341 »  340 ,  341 

340  1F(IS  ,ttT.  1)  60  TO  321 

341  CONTINUE 

IF (N  “D  302,301,302 

301  EAH (1*1) =£Ap  ( A ( 1 f 1 ) *H ) 

£2 Ah (1*1) =c  AH  (1*1) *E AH (1,1) 

E3Ah  (1,1)  -E2AH  (1*1)  *EAH  (1,1) 

AH2  ( 1 » 1 ) -AH (1*11  * AH (111) 

IF( INDEX  ,lu.  1)  60  TO  350 
AH3<1»1)=1./(AH2(1,1)*AH(1»1>  ) 

60  TO  303 

302  CONTINUE 

DO  330  I— 1 , N 
DO  331  J=1 ,  N 
DO  332  K=1»N 

AH2(I,J)=AH2(I»J)+AH(I»K)*Ah(K,J> 

332  CONTINUE 
331  CONTINUE 
330  CONTINUE 

DO  333  I=1*N 
DO  334  J=1*N 
DO  335  K=1*N 

EAH( I » J) =EAH (  I ,  J) +AH2  (I,K)*AH(K,J) 

335  CONTINUE 
334  CONTINUE 

333  CONTINUE 

IF ( INDEX  .EU.  1)  60  TO  351 
CALL  INVERT(EAH,N,AH3) 
call  PADE(A»H*EAH,N> 

H2=H*H 

CALL  PADE(A,H2»E2AH,N) 

H3=H2+H 

CALL  PACE ( A , H3 , E3AH » N ) 

303  CONTINUE 

DO  304  1=1 ,N 
DO  305  J=1»N 

M  ( 1  *  1 » I »U)=UnIT ( I # J) -1 «5*AH( I » J)+AH2 ( I  #  J) 

W  ( 1  •  2 » I » J)=UNIT ( I , J) -0«5*AH( I ,  J) 

*<1.3,I.J)=UNIT(I,J)+0.5*AH(I,J> 

ft  ( 1  •  4*  I  *  J)=UNIT  ( 1 , J) +1 «5*AH( I , J)+AH2(  I,  J) 

W  ( 2«  1 , 1 1  J)  =-2  •  *  ( UNIT  ( I  *  J) -Ah(  I » J) ) 

W(2,2»I,JJ=-2.*UnIT(I, J)+AH2(I*J) 
tf (2 t 3,  I » J)=-2«* (UNIT ( I , J) +AH( I » J) > 
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W<2*4*  I#U)=-E.*UNIT(I#J)-4.*AH(I»J)-3.*AH2(I#  J) 

W(3,1»I#J)=UNITU»  J)-O.S*AH(I,J) 

M  (3f  2 » I » J)=L'NIT  ( I » J)  +0.5* AH  ( I » J) 

W(3,3.I.J)=Ml#4#I,J) 

W(3,4#I*U)=UNIT(I*  JM-2.5*AH(I,J)+3,*AH2(I.J) 

305  CONTINUE 
304  CONTINUE 

360  DO  306  U=i#KUP 
00  307  I=i»N 
00  308  J=1>N 
00  309  K=1#N 

QHI  <  1 1  ■'  I  r  J)  =tlHI  1 1 1 » I » J)  ALPHA  (1)*W(II»1*I  *K)  *E3A*~*  ( K  #  J) ♦ ALPHA  (2 ) *W ( 
SI  I  *2»  I  *  K )  *c.EAH(K#  J) +ALPHA  { 3 )  ( 1 1  *3»  I  *K  )*EAH(K » J) 

309  CONTINUE 

QHl(Il»I#J)=QHICII,I,J)*W<Ii.4.I*J> 

308  CONTINUE 
307  CONTlNUt 

306  CONTlNUt 

00  310  K=i'KUP 
00  311  1=1»N 
00  312  J=1'N 
00  313  L=1>N 

IF(INOEX.EO.O)  PHI  'K»  I  *  J)sPHI  (K»  I » J)-AH3(  I»L)*QhI  (K»L»  J) 

IF ( INDEX. Eu.l)  PH:(KfI»J)=PHI(K»l»J)-AH4(I*L)*QHl(K,L,J) 

313  CONTINUE 
312  CONTINUE 
311  CONTINUE 

310  CONTINUE 
321  CONTINUE 

UO  314  K=1»KUP 

CALL  6FN(6#MiN»Y»K#T#A) 

00  315  I=1*N 
00  316  J=1 t N 

YN( I ) =YN( I )*H*PHI  (KtI*J)*G(j) 

IF ( K  »EQ«  KUP)  YN(I)=YN(I)-ALPHA(3)*EAH(I»J)*Y(3*d)-ALPHA(2)*E2AH( 
SI » J) *Y (2»J) “ALPHA  ( 1 ) *E3AH( I  * J) *Y ( 1 » J) 

316  CONTINUE 
31b  CONTINUE 

314  CONTINUE 
RETURN 

350  AH3 (1*1) = AH2 ( 1  r  1 )  *  AH  ( 1 » 1 ) 

AH4 ( 1 > 1 ) =1 <  0/ ( AH3  (1»1) *AH (1*1)) 

60  TO  357 

351  00  352  IsltN 
00  353  J=1*N 
AH3( I t J)SEAH(ItJ) 

EAH(I»J)=0.0 

353  CONTINUE 

352  CONTINUE 

00  354  1=1 rN 
00  355  J=liN 
00  356  K=1*N 

EAH(  I  #  J)=EAh(  I#U)+AH3(  I  »K)  *AH(K»  J) 

356  CONTINUE 
355  CONTINUE 

354  CONTINUE 

CALL  INVERT (EAH»N»AH4| 

CALL  PACE (AiH*EAH»N) 
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357 


359 

358 


H2=H+H 


CALL  PA0E(A#H2#t2AH#N> 
H3=rt2*H 

CALL  PAD£ ( A » H3 » L 3 AH #  N ) 
CONTINUE 
UO  358  I=1»M 
UO  359  J=1»N 


WU 
w  ( 1 
mi 
mi 
¥  {?. 
w(2 
In  (2 
W<2 
W(3 
w(  3 
M  (3 
W(3 
rf(4 
W  ( 4 
to  ( 4 
W  ( 4 


.J) 

rJ) 

#U) 

#J> 

,J) 

#J> 

fJ) 

,J) 

,J) 


CONTINUE 

continue 

60  TO  360 
END 


-0.1*UN1T< 
-0.1*UNIT( 
-0.1*UNIT( 
-0,l*UNIT< 
0.3*UNIT(I 
0<3*UnIT( I 
0.3*UNIT(I 
0.3*UNIT(I 
-0,3*UN1T( 
-0,3*UNIT( 
-0t3*UNIT( 
-0,3*UNIT< 
O.lMUNITt 
0.1*UNIT(I 
0« 1* (UNIT ( 
0.1*UNIT(I 


I#J)4l.l*AH(I#J)-23.*AH2(I,J)/15.4AH3(I»d> 

I»J)*AH(I»d)-29,*AH2<I»J>/60. 

I  * J) ♦0.9*AH ( I  »U)+7,*AH2 ( I r  J) /15. 
IiJ)*0.6*AH(I'J)*79.*AH2(I'J)/60.  +  .9*AH3(I'J) 
»J)-2t3*AH(I»J)*2.1*AH2(I» J) 
»J)-2.*AH(I»J)-0.05*AH2(If J)+AH3(I.J) 
'J)-1.7*AH(I'J)-1.9*AH2(I'J> 
#J)-1.4*AH(I»J)-3.45*AM2(I,J)-2.7*AH3(I»J) 

I » J) +1 . 3*AH ( I#  J)-0.6*AH2(  I » J) 

I»J)*Ah(X'J)  +  .55*AH2(I’J) 

IfJ)+0,7*AH( I r J) +1.4*AH2  ( I #  J)+AH3( I » J) 
I»J)>0.4*AH(I»J)+1.95*AH2(I»J)+2.7*AH3(I*J) 

X»  J)-AH(I»J))4AH2(I»J)/30. 

.J)-AH2(I* J)/60, 

I*J)4-Ah(I*J)  )*AH2(I» Jl/30. 

»J)+0«2*AH(IrJ)  ♦11,*AH2  ( I » j)  /60  .+0. 1*AH3  ( I » J) 


LISTS  OF  SYMBOLS  AND  DEFINITIONS 


Indices 

i,  j,  K,  0. ,  m,  n,  N,  p,  q,  v 

Scalars 

G,  E,  k,  K#,  L,  L*,  L#,  Q#,  3 
a  ,  a  ,  0  ,  Y  ,  5  (h) ,  <  ,  ,  0  ,  <  ,  £  ,  *  ,  4>  ,  T ,  A 

Vectors 

e  ,  f,  g  ,  G ,  w ,  X  y  ,  y',  y*,  *,4,,^,  eg 

Matrices 

A,  C(h),E,H,K,Q,  4,  ^Ki (Ah) 

Symbols 
! 

C 

I  I 

II  II 

» 

3 

? 

V 

n 

£ 

LHS.RHS 


factorial 
belongs  to 
absolute  value 
norm 

much  greater  than 
there  exists 
such  that 
for  every 
product 
summation 

left-hand  side,  right-hand  side 
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'V 
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LMS 

NLMS 

t 

t' 

h 

hN 

i 

PN 


V«n’  *n 


0) 


9 


#.  * 

p  (t) 


(Pp 

G(y) 

X]  (Ah) 

1 

i 

gj(0 

X(A) 

p(A) 

p(f).p(D.p(x,n 

0(hP+2) 

„p+l 


Linear  Multistep 
Nonlinear  Multistep 
time  variable 
integration  variable 
step  size 

step  size  of  method  N 

roots  of  characteristic  polynomial 

method  N  of  order  p 

initial  vector 

values  of  y ,  g ,  f  evaluated  at  t 

j-th  derivative  of  g 

a  polynomial  in  r 

class  of  polynomials  of  degree  p 

a  function  of  y 

an  integral 

a  function  o  flj  (Ah) 

M.  O.  C.  (modulus  of  continuity) 
eigenvalues  of  A 
spectral  radius  of  A 
characteristic  polynomials 
order  of  hP+2 

functions  which  have  (p+l)-th  continuous  derivative 


error  function 


cCow] 

r[y(t);h] 


true  operator 
nonlinear  operator 


local  discretization  error 
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